^■ad-ami  551  INVESTIGATION  OF  THE  MECHANISM  Of  TRANSONIC  SHOCK  1// 

WAVE/BOUNDARY  LAYER  INT..(U)  SCIENTIFIC  RESEARCH 
ASSOCIATES  INC  GLASTONBURY  CT  D  V  ROSCOE  ET  AL .  MAY  84 
UNCLASSIFIED  SRA-R83'930006-F  ARO- 1 7048 . 1 -EG  F/G  20/4 


NL 


MICROCOPY  RESOLUTION  TEST  CHART 
national  eur<tAU  of  SrANDAROS19&3  a 


Ajlo  non.t*^ 


AD-A141  551 


INVESTIGATION  OF  THE  MECHANISM  OF  TBANSONIC  SHDCE  NAVE/BODNDAKY 
LATER  INTERACTIONS  USING  A  MAVIER-ST(»ZS  ANALYSIS 

(FINAL  REPORTl 


D.V*  Roscoc 
S.J.  Shaoroth 
H.J.  Glbellng 
H.  McDonald 


May  1984 


U.S.  ARMY  RESEARCH  OFFICE 
Contract  DAA629-80-(M)082 


Freparcd  by 

SCIENTIFIC  RESEARCH  ASSOCUTES,  INC. 
P.O.  Box  498 
Glastoobnty*  Cl  08033 


C 


DTIC 

ELECTE] 

MAYaSHM 


SeeURlTV  CLASSirtCATlON  OF  THIS  PAGE  (Whmt  0«r«  Snfemd) 


•7" 


•  i 


REPORT  DOCUMENTATION  PAGE 

READ  mSTRUCnONS 

BEFORE  COMPLETINC  FORM 

1.  ni^ohf  NUMOER  2.  OOVT  ACCESSION  NO. 

N/A 

3.  REClPtCNT'S  CATALOG  NUMBER 

N/^ 

4.  title  ranN  SuMtlaj 

Investigation  of  the  Mechanism  of  Transonic 
Shockwave/Boundary  Layer  Interactions  Using 
a  Navier-Stokes  Analysis 

S.  TYRE  OR  REPORT  h  PERIOD  COVERED 

Final  Report 

July  1980  -i.tor7l983 

«.  PERPORMINO  ORO.  REPORT  NUMEER 

SRA  Report  R83-930006-F 

7.  AUTNOR(a> 

D.  V.  Roscoe,  S.  J.  Shamroth, 

H.  J.  Gibeling,  H.  McDonald 

S.  CONTRACT  OR  grant  NUMEERfaf 

DAAG29-80-C-0082 

S.  reilPORMINO  OROANIZATION  NAME  AND  ADDRESS 

Scientific  Research  Associates,  Inc. 

P.O.  Box  A98 

100  Sycamore  Street 

Glastonbury.  CT  06033 

to.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  0  WORK  UNIT  NUMBERS 

II.  CONTROLLINO  or  PICE  NAME  AND  ADDRESS 

U.  S.  Army  Research  Office 

Post  Office  Box  12211 

Researrh  Tr-lanolo  Pai-V  Mr  OTTno 

12.  REPORT  DATE 

May  1984 

I*.  NUMBER  OP  PAGES 

14.  MONITORINO  AOENCy'fiAME  4  ADDRESSCIf  Mllmniti  Mai  CanMItin*  Olllea) 

IS.  SECURITY  CLASS,  (oi  Aim  rwpott) 

Unclassified 

18a.  OECLASSIFICATION/OOWNGRAOINC 

schedule 

1#.  OlSTfnOUTION  STATCMCNT  (•t  thim  AMporl} 

Approved  for  public  release;  distribution  unlimited. 

17.  OISTRISUTION  statement  (•!  (ha  ahalraet  aalara.  In  Sloe*  20.  If  dlllmnni  Mai  Rmpowt) 

NA 

IS.  surruememtary  notes 

The  view,  opinions,  and/or  findings  contained  in  this  report  are 
those  of  the  author(s)  and  should  not  be  construed  as  an  official 

Department  of  the  Army  position,  policy,  or  decision,  unless  so 

10.  KEy  WOROS  (CnnUnim  aa  nwm—  alM  II  aacaaaaiT  mtd  Idtntllr  hr  Maeh  maaharf 

Multidimensional  Implicit,  Time  Marching,  Navier-Stokes,  Shock  Wave/Boundary 
Layer  Interaction,  Adaptive  Grid,  Transonic  Bump  Flow,  Impinging  Shock, 

Normal  Shock 

.  1. 

2E.yaeET1IACrrCha*MM  anavaraa  m  aaaaaaair  mE  MmHllr  hr  hlaah  wharf 

The  present  report  discusses  the  development  of  a  time-dependent  Navier- 
Stokes  code  for  use  in  predicting  transonic  shock-wave  boundary  layer  interaction 
In  addition,  various  test  cases  which  have  been  performed  are  discussed.  The 
algorithm  used  to  solve  the  equations  is  based  upon  the  consistently  split 
linearized  block  implicit  method  of  Briley  and  McDonald  [15,  16].  The  philosophy 
and  use  of  a  solution  adaptive  mesh  is  also  described.  The  test  cases  studied 
in  this  report  are:  normal  shock  wave  boundary  layer  interaction  in  a  constant  — 

D0,;2r»Un 

UNCLASSIFIED 


■Kcumrr  ei.«UiriCATION  or  TMIS  r*«K  (rnmn  Oata  Batara« 
1 


f 


UNCLASSIFIED 


ttCURITV  CUMUFICATIOH  OP  THIS  PMtfWIiwi  Oat*  BntwaO 


11 


LIST  OF  SYMBOLS 


van  Driest  damping  coefficient 
specific  heat  at  constant  pressure 
determinant  of  the  Jacobian  matrix 
dissipation  function 
distance  to  the  nearest  wall 
dimensionless  distance  to  the  nearest  wall 
enthalpy,  throat  height 
mixing  length 

mixing  length  in  the  core  flow  region 
static  pressure 

magnitude  of  the  velocity  / 


turbulent  heat  flux  vector 

mean  heat  flux  vector 
universal  gas  constant 
Reynolds  number 
time 

temperature 

stagnation  temperature 

velocity  vector 
velocity  component  in  x-direction 
friction  velocity 
velocity  component  in  y-direction 


Accession  For 

”nt:S  GRA4I 

DTIC  TAB 

□ 

Unamioanced 

n 

Justification — 

■  ■■ 

i 


Hv  _  _ 

DlaW 

Aval 

Dlst 

.butloa 

lablllt 

Avail 

Spac 

/  _ _ i 

y  Codes 

and/or 

lal 

1 

111 


LIST  OF  SYMBOLS  (Continued) 


Symbols 

w  velocity  component  in  z-direction 

Wg  w  at  the  edge  of  the  bomidary  layer 

X,  XX  cartesian  coordinate  in  transverse  direction 

y,  y2  cartesian  coordinate  in  spanwise  direction 

y^  computational  coordinates 

z,  X3  cartesian  coordinate  in  streamwise  direction 


Greek  Symbols 
6 


e 

K 


V 

''art 

5,  n,  ? 


boundary  layer  thickness 
turbulence  energy  dissipation  rate 
von  Karman  constant 
dynamic  viscosity 
artificial  dissipation 
computational  coordinates 


X 

xT 

P 

o 

T 

"^xx*  '^xy» 

♦ 


molecular  stress  tensor 
turbulent  stress  tensor 
density 

artificial  dissipation  parameter 
time 

local  shear  stress 
component  of  stress  tensor 
meanflow  dissipation  rate 


Iv 


LIST  OF  SYMBOLS  (Continued) 

associated  with  the  bottom  wall 
associated  with  the  side  wall 
associated  with  the  time  or  top  wall 
associated  with  the  x-direction 
associated  with  the  y-direction 
associated  with  the  z-direction 

associated  with  turbulent  quantities, 
transpose  of  matrix 


TABLE  OF  CONTENTS 


INTRODUCTION  . 

ANALYSIS  . 

Governing  Equations  . 

Turbulence  Modeling  . 

Zero  Equation  Model  -  (Mixing  Length)  . 

One-Equation  Model  -  (k-4)  . 

Two-Equation  Model  -  (k-c)  . 

COORDINATE  SYSTEM  DETAILS  . 

Solution  Procedure  . 

DISCUSSION  OF  RESULTS  . 

Normal  Shock  Wave-Boundary  Layer  Interaction  in  a 

Constant  Area  Tube  . 

Boundary  Conditions  . 

Solution  Adaptive  Coordinate  System  .  .  .  . 
Results  . 

Steady  Transonic  Flow  Over  an  Axisymmetric  Bump  . 

Boundary  Conditions  . 

(i)  Large  Radius  of  Curvature  Bump; 

Thin  Inlet  Boundary  Layer  . 

(ii)  Large  Radius  of  Curvature  Bump; 

Thick  Inlet  Boundary  Layer  . 

(iii)  Small  Radius  of  Curvature  Bump; 

Thin  Inlet  Boundary  Layer  . 

Steady  Transonic  Flow  Over  an  Axisymmetric 
Bump:  Sunnary  . 

Unsteady  Transonic  Flow  Over  an  Axisymmetric  Bump 
(i)  High  Frequency  Case  (fj.  ■  0.175).  .  .  . 
(ii)  Low  Frequency  Case  (f^  *  0.004)  .  .  .  . 

Unsteady  Flow  Cases:  Suimnary  . 

Steady  Flow  Over  an  Axisymmetric  Compresion  Ramp 

Oblique  Shock-Wave  Impingement  on  a  Flat  Plate 
Boundary  Layer  . 

REFERENCES  . 

APPENDIX  -  SOLUTION  PROCEDURE  [17]  . 

Background  . 

Split  LBI  Algorithm  . 

Linearization  and  Time  Differencing  .  .  .  . 
Special  Treatment  of  Diffusive  Terms  .  .  .  . 
Consistent  Splitting  of  the  LBI  Scheme  .  .  . 


FIGURES 


1 


INTRODUCTION 


The  flow  field  resulting  from  the  interaction  of  a  shock  wave  and  a 
boundary  layer  on  a  projectile  or  missile  remains  a  major  problem  which  has 
yet  to  be  completely  analyzed.  The  interaction  process  plays  an  important 
role  in  several  areas.  Among  these  are:  (i)  body  lift  and  drag,  (ii)  side 
force  development  and  control  and  (iii)  heat  transfer.  Since  the  interaction 
causes  an  abrupt  pressure  rise  and  boundary  layer  thickening  and  may  be  ac¬ 
companied  by  regions  of  local  separation,  the  interaction  flow  can  be  a  major 
contributor  to  the  overall  body  drag  and  can  cause  substantial  changes  from 
the  potential  flow  pressure  distribution.  Furthermore,  since  the  interac¬ 
tions  are  often  three-dimensional  due  to  either  geometric  or  flow  incidence 
effects,  the  generated  forces  may  not  be  symmetric  and  may  result  in  sigifi- 
cant  side  forces.  These  side  forces,  in  turn,  can  lead  to  serious  control 
and  guidance  problems.  Finally,  the  shock  wave  boundary  layer  interaction 
zone  may  be  a  region  of  severe  heat  transfer.  The  problems  associated  with 
interactions  are  particularly  troublesome  at  transonic  speeds  where  both  the 
shock  location  and  its  shape  are  very  sensitive  to  minor  changes  in  flow 
geometry.  As  a  consequence,  in  recent  years  there  has  been  considerable 
interest  in  the  development  of  accurate  and  efficient  prediction  techniques 
for  this  category  of  the  interaction  problem. 

At  present,  two  main  approaches  for  treating  the  shock  wave-boundary 
layer  interaction  problem  are  being  pursued.  The  first  termed  the  strong  in¬ 
teraction  approach  divides  the  flow  region  of  interest  into  two  parts,  an 
outer  inviscid  part  and  a  wall  viscous  part.  Each  region  is  then  solved  sep¬ 
arately  via  the  appropriate  set  of  equations  for  that  region.  Equations  ap¬ 
propriate  for  inviscid  flows;  e.g.,  Euler,  potential  flow  or  simple  wave  re¬ 
lations,  are  used  in  the  outer  region  and  boundary  layer  equations  are  used 
in  the  inner  region.  At  the  juncture  of  the  regions,  matching  conditions 
which  require  continuity  of  all  flow  variables  are  applied.  In  flow  situa¬ 
tions  in  «diich  the  outer  flow  is  supersonic  and  thus  described  by  hyperbolic 
equations,  a  solution  can  be  forward-marched  in  space  with  the  inviscid  and 
viscous  regions  coupled  on  a  station-by-station  basis.  The  chief  difficulty 
in  this  process  is  that  mathematically  stiff  equations  must  be  solved.  Com¬ 
mon  problems  with  stiff  equations  manifest  themselves  in  the  form  of  numeri¬ 
cal  solutions  which  can  branch  off  the  desired  solution  thus  producing  a 
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physically  unrealistic  result.  In  regions  where  the  inviscid  flow  is 
subsonic  and  thus  described  by  elliptic  equations,  a  forward-marching 
procedure  without  iteration  is  not  physically  realistic  since  the  upstream 
pressure  propagation  is  not  modeled.  Consequently  a  sequence  of  inviscid  and 
boundary  layer  solutions  must  be  performed  in  a  manner  in  which  each  stage 
corrects  the  former  stage  through  global  iteration.  Additional  problems 
occur  in  transonic  flows  where  both  supersonic  and  subsonic  outer  regions  are 
present  and  where  small  displacement  effects  may  have  considerable  influence 
on  shock  location  and  overall  pressure  distribution.  In  cases  where  the 
nominally  inviscid  flow  is  subsonic  behind  the  shock  the  situation  is  further 
complicated  by  the  subsonic  outer  region  being  elliptic.  Since  the  flow 
cannot  be  forward  marched  here,  a  global  rather  than  a  station-by-station 
iteration  must  be  used. 

If  an  iteration  procedure  is  to  be  used,  the  viscous  layer  can  be 
solved  iteratively  by  either  a  forward-marching  boundary  layer  calculation 
procedure  or  as  the  asymptotic  condition  of  a  time  or  pseudo-time  dependent 
integration.  In  the  case  of  steady  state  forward-marching  boundary  layer 
procedures,  problems  are  encountered  when  the  boundary  layer  is  subjected  to 
an  adverse  pressure  gradient  strong  enough  to  cause  separation.  Although  a 
boundary  layer  procedure  can  be  forward  marched  through  separation  by 
suppressing  streamwise  convection  terms  in  the  separated  flow  region,  via 
the  so  called  FLARE  approximation  [1],  the  resulting  solution  is  based  upon 
an  approximation  made  in  the  separated  flow  region  which  becomes 
progressively  more  inaccurate  as  the  backflow  velocities  increase. 

Therefore,  calculated  details  of  the  flow  in  this  region  cannot  be  expected 
to  be  accurate  with  significant  backflow  using  FLARE.  A  global,  but  time 
consuming,  iteration  [2]  is  necessary  to  replace  the  FLARE  approximation,  and 
at  this  point  the  efficiency  of  the  forward-marching  scheme  must  be  carefully 
evaluated  to  ensure  a  net  gain  exists  relative  to  solving  the  full 
Navier-Stokes  equations  for  the  viscous  layer.  Time  integration  of  the 
interaction  between  the  boundary  layer  and  the  external  flow  can  be 
structured  to  avoid  the  use  of  either  FLARE  or  the  global  iteration  to 
account  specifically  for  the  backflow  velocity  [2],  and  in  some  cases  the 
time  marched  iteration  approach  may  show  a  significant  gain  relative  to 
solving  the  Navier-Stokes  equations  for  the  viscous  flow  region.  However, 
the  interaction  approach  remains  limited,  even  with  time  integration,  to 
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flows  where  the  division  of  the  flow  into  zones  which  can  be  interacted  is 
reasonable  and  details  of  any  separated  zone  are  not  of  major  interest.  A 
discussion  of  some  of  the  early  techniques  aimed  at  the  shock  wave-boundary 
layer  interaction  problem  which  are  based  upon  a  viscous-inviscid  approach  is 
given  by  Hankey  [3] . 

The  second  approach  used  in  obtaining  solutions  for  shock  wave-boundary 
layer  flow  field  is  the  fully  viscous  analysis.  The  fully  viscous  analysis 
solves  the  entire  flow  field  via  a  single  set  of  equations,  thus  avoiding  any 
need  to  divide  the  flow  into  viscous  and  inviscid  regions.  The  Navier-Stokes 
analyses  pursued  to  date  include  explicit,  implicit  and  hybrid  numerical  for¬ 
mulations.  For  example,  in  a  very  early  work  of  this  type,  MacCormack  [4] 
applied  an  explicit  procedure  to  the  laminar  interaction  problem.  In  subse¬ 
quent  work  MacCormack  and  his  co-workers  have  used  this  explicit  method 
extensively  in  analyzing  the  interaction  flow  field  [5-8).  The  basic  method 
has  also  been  used  by  Shang  and  Hankey  [9]  in  calculating  turbulent  flow  over 
a  compression  ramp  and  in  predicting  a  streamwise  corner  interaction  region 
[10].  Recently  MacCormack  has  developed  a  new  hybrid  explicit-implicit 
characteristics  method  [11]  which  is  considerably  faster  than  the  fully 
explicit  method  [4-10],  and  has  applied  this  technique  to  shock  wave-boundary 
layer  interaction  problems.  This  same  technique  has  also  been  applied  to 
axisymmetr ic  interaction  flow  fields  by  Viegas  and  Coakley  [12]. 

In  addition  to  the  fully  explicit  investigations  of  Ref.  4-10,  and  the 
hybrid  investigations  of  Refs.  11  and  12,  Levy,  Shamroth,  Gibeling  and 
McDonald  [13]  and  Beam  and  Warming  [14]  have  applied  implicit  solution 
procedures  to  the  interaction  flow  field  problem.  In  the  former  invest¬ 
igation,  Levy  et  al  applied  the  Briley-McDonald  consistently  split  Linearized 
Block  Implicit  (LBI)  solution  procedure  [15,  16]  to  the  shock  wave-boundary 
layer  interaction  problem  in  a  feasibility  study  which  was  part  of  a  larger 
effort  primarily  focused  on  a  boundary  layer  strong  interaction  study.  Beam 
and  Warming  [14]  also  applied  a  similar  implicit  scheme  to  the  interaction 
problem. 

To  date  many  investigations  have  focused  upon  the  interaction  flow 
field  problem.  Most  of  the  analyses  have  concentrated  upon  two  specific 
problems;  ( i)  flat  plate  -  ramp  compression  corners  and  (ii)  two-dimensional 
shock  impingement  on  a  flat  plate.  The  problem  of  axisymmetr ic  transonic 
interactions  which  may  be  particularly  relevant  in  missile  or  projectile 
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applicaCions  has  received  less  attention.  Prior  to  the  present  study,  Viegas 
and  Coakley  [12]  applied  MacCormack's  hybrid  technique  to  an  interaction 
occurring  on  an  axisymmetr ic  body.  Little  other  work  has  concentrated  upon 
the  axisymnetric  problem. 

This  report  presents  and  describes  the  development  and  application  of  a 
fully  viscous  method  for  solving  transonic  shock  wave-boundary  layer 
interaction  problems.  The  method  utilizes  the  Briley  and  McDonald  1^13,  16] 
linearized  block  implicit  (LBI)  technique  for  the  numerical  solu'  .i  jf  the 
multidimensional,  time-dependent  Navier-Stokes  equations,  and  ha  ?en 
applied  to  several  axisymmetr ic  flow  problems.  These  include  tr  onic  flow 
in  a  constant  area  axisymmetric  pipe,  transonic  flow  over  an  axi  jia  .ric 
bump,  shock  inducement  at  an  axisyrmnetric  compression  corner  and  tne 
two-dimensional  flat  plate  boundary  layer  -  impinging  shock  wave  interaction 
problem.  In  addition  to  these  steady  flow  cases,  several  time-dependent 
problems  are  reported  in  which  the  transonic  flow  over  an  axisymmetric  bump 
was  subjected  to  downstrean  static  pressure  oscillation. 

In  the  following  sections,  the  governing  equations  are  presented  and  a 
general  description  of  their  transformation  to  general  curvilinear 
coordinates  is  given.  However,  description  of  the  specific  coordinate 
transformation  used  for  each  case  is  presented  on  a  case-by-case  basis. 


ANALYSIS 

Governing  Equations 

The  equations  used  in  the  present  effort  are  the  ensemble-averaged, 
time-dependent  Navier-Stokes  equations  which  can  be  written  in  vector  form  as 
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Momentum 
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Energy 


+  V  •  (puh)  =  -  V  •  (q  -t-  q^)  +  ^  +  *  +  pe  (3) 

where  p  is  density,  u  is  velocity,  p  is  pressure,  n  is  the  molecular  stress 
tensor  is  the  turbulent  stress  tensor,  h  is  enthalpy,  q  is  the  mean  heat 
flux  vector,  q"^  is  the  turbulent  heat  flux  vector,  ♦  is  the  mean  flow 
dissipation  rate  and  e  is  the  turbulence  energy  dissipation  rate.  If  the 
flow  is  assumed  at  a  constant  total  temperature,  the  energy  equation  is 
replaced  by 

2 

^  “  constant  (4) 

o  2C 

P 

where  Tq  is  the  stagnation  temperature,  q  is  the  magnitude  of  the  velocity 
and  Cp  is  the  specific  heat  at  constant  pressure.  In  a  number  of  cases 
considered  in  this  work  the  assumption  of  constant  total  temperature  has  been 
invoked  by  using  Eq.  4  as  an  approximation  to  Eq.  3  for  the  sole  purpose  of 
reducing  computer  run  time  where  the  constant  Tq  assumption  was  warranted. 
Cases  in  which  this  substitution  has  been  made  are  identified  in  the  text 
describing  results.  A  number  of  terms  appearing  in  Eqs.  1-4  require 
definition.  The  stress  tensor  appearing  in  Eq.  2  is  defined  as 

It  -  2pD  -  (-|  u-Kg)  7  •  C  n  (5) 

where  Kg  is  the  bulk  viscosity  coefficient  and  is  the  deformation  tensor, 
defined  by; 

D  =  Y  ((VU)  +  (VU)'^)  (6) 

In  addition  the  turbulent  stress  tensor  has  been  modeled  using  an  isotropic 
eddy  viscosity  such  that: 


-  p  u'  u'  -  2p.j,ID  (7) 

where  ii-f*  the  turbulent  viscosity,  is  determined  by  a  suitable  turbulence 
model.  Turbulence  modeling  is  described  in  some  detail  in  the  next  section. 
Equation  3  contains  a  mean  heat  flux  vector  defined  as  follows: 


and  a  turbulent  heat  flux  vector  defined  as: 

q'*’  -  (7t) 

where  <  and  are  the  mean  and  turbulent  thermal  conductivities, 
respectively. 

Also  appearing  in  Eq.  3  is  the  mean  flow  dissipation  term 
«  -  2uID  :  id  -  (4  W  -  K^)  (V  •  U)^ 


The  equation  of  state  for  a  perfect  gas 

p  »  pRT 

where  R  is  the  gas  constant,  the  caloric  equation  of  state 

e  ■  c  T 

V 

and  the  definition  of  static  enthalpy 

h  -  c  T 
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(9) 


(10) 


(11) 


(12) 


(13) 


supplement  the  equations  of  motion. 

Finally  the  flow  properties  u,  k  and  Rg  are  determined  using  the  following 
constitutive  relations. 

The  molecular  viscosity  p  is  determined  using  Sutherland's  law. 

^  3/2  T  +  S. 

^  •  (I-  )  ‘ 


V  '•T 
o  o 


T  +  S 


(14) 
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where  Sj  ■  100*K  for  air. 

The  bulk  viscosity  will  be  assumed  to  be  zero. 


(15) 


and  the  thermal  conductivity  may  be  determined  by  use  of  a  relation  similar 
to  Sutherland's  law  viz. 


^  _  3/2  T^  S, 

—  m  fZ_  )  O  Z 

K  Ij  J  T  ♦  S« 
o  o  2 


(16) 


where  S2  ■  194*K  for  air. 
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Turb;.  ?nce  Modeling 


Several  alternative  turbulence  models  have  been  considered  in  the  | 

course  of  the  present  work.  In  general  terms,  the  models  used  were  zero,  f 

one-  and  two-equation  models.  The  formulation  of  each  is  described  in  this 
section . 

Zero  Equation  Model  -  (Mixing  Length) 

Of  all  available  turbulence  models,  Prandtl's  mixing  length  model  is 
probably  still  the  most  widely  used.  The  model  was  originally  developed  for 
use  in  unseparated  boundary  layer  flow  situations  and  has  been  shown  to 
perform  well  under  such  conditions.  In  the  cases  described  in  this  report 
the  mixing  length  model  has  been  used  extensively  in  order  to  investigate  its 
use  in  the  shock  wave/boundary  layer  interaction  problem.  An  advantage  of 
the  method  from  the  point  of  view  of  economy  is  that  is  does  not  require 
additional  transport  equations  to  model  the  effect  of  turbulence,  but  rather 
relates  the  Reynold's  shear  stress  to  mean  flow  quantities  via: 

PuTuT  ■  ^T 

"  J  R- 
e 


where 


where 


1J.J.  -  (2D  :  D 

I  ■  min[i  ,  <dD] 


where  d  is  the  normal  distance  to  the  nearest  wall  and  D  is  van  Driest 
damping  coefficient  given  by 

D  ■  1  -  exp(-  y  /A  ) 

-  0.096 

<  -  0.4 
♦ 


(17) 


y  ■  du^/v 


1/2 


where  the  local  shear  stress  tji  is  obtained  from 

>1/2 


-  (2D  :  D  )' 


(18) 


and  D  is  defined  by  Eq.  6. 


One  problem  in  Che  mixing  length  formulation  is  Che  definition  of  5.  In 
boundary  layers  the  sCreamwise  velocity  u  approaches  an  edge  velocity,  Ug, 
asymptotically,  however,  a  monoConic  approach  to  an  asympoCoCic  edge  velocity 
is  not  characteristic  of  Navier-Stokes  solutions.  In  order  to  avoid  the 
problem  of  determining  Che  boundary  layer  edge,  5,  as  defined  in  Che  usual 
boundary  layer  context,  i.e.,  6  is  the  distance  from  Che  wall  at  which  u/ug 
>  0.99,  Che  following  relation  is  used. 
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2.0d 


(9/9max”c) 


(19) 


In  other  words,  d,  is  taken  as  twice  Che  distance  (measured  away  from  Che 
nearest  wall)  for  which  q/q^ax  “  *-•  value  of  c  used  in  the  present 

effort  was  0.90. 


One-Equation  Model  -  (k-t) 

Although  as  discussed  above  Che  mixing  length  concept  is  valid  for  a 
variety  of  flows,  some  important  flow  situations  arise  in  which  a  less 
restrictive  model  is  required.  One  such  model  is  the  one-equation  turbulence 
model  (17,  18,  19]  in  which  a  transport  equation  for  turbulence  kinetic 
energy,  k,  is  formulated  as  follows: 

Iy—  ♦  V»(pUk)  -  7k)  +  Iv-i  (  D:  D)  -  pe 

k 

where  following  Ref.  20,  Oj^  is  set  to  1.0,  and  k  is  the  turbulence  kinetic 

energy  _ 

k  »  Y  u'  •  u^  (20) 

and  Che  Pr and cl -Kolmogorov  relaCion,  Eq.  21 ,  defines  Che  CurbulenC  viscosiCy 

as: 


In  addition,  Che  turbulence  dissipation  rate  e  is  determined  by: 

c . 


(21) 


(22) 


where  A,  is  a  relevant  turbulent  length  scale  for  the  problem  of  interest. 
The  k-t  model  has  an  advantage  over  the  mixing  length  model  in  Chat  Che  use 


of  a  transport  equation  for  turbulence  kinetic  energy  allows  for  its 
convection,  production  and  dissipation.  This  is  important  because  it  allows 
a  nonequilibrium  effect  on  the  turbulence  to  be  included  in  the  calculation 
while  the  mixing  length  model  can  only  account  for  local  equilibrium 
turbulence  effects  via  its  association  with  the  mean  velocity  field.  A  major 
disadvantage  which  the  k-i  model  shares  with  the  mixing  length  model  is  the 
requirement  of  length  scale  specification.  Typically  the  mixing  length,  as 
described  in  the  previous  section,  is  used  as  the  representative  length 
scale . 

In  modeling  the  flow  in  the  near  wall  region  where  low  local  turbulence 
Reynolds'  numbers  occur,  two  approaches  are  available.  The  first  is  the  wall 
function  approach  which  does  not  resolve  the  near  wall  region  but  assumes 
specific  functional  forms  for  the  required  turbulence  quantities  and  uses 
these  forms  to  create  the  required  normal  derivative  formulations  at  the 
first  grid  point  from  the  wall.  Such  an  approach  obviously  requires  a 
detailed  knowledge  of  the  turbulence  model  dependent  variables  in  the 
vicinity  of  the  wall.  Although  reasonable  functional  formulations  can  be 
specified  for  simple  flows  such  as  constant  pressure  boundary  layers, 
specification  in  the  much  more  complex  flows  of  current  interest  are  much 
more  difficult.  Therefore,  the  alternative  approach  in  which  the  viscous 
sublayer  is  resolved  has  been  used.  The  method  makes  no  approximation  at  the 
boundary,  but  requires  that  the  near  wall  low  turbulence  Reynolds'  number 
physics  be  modeled.  In  this  work,  a  near  wall  model  which  was  successfully 
used  by  Shamroth  and  Gibeling  [21]  in  a  time-dependent  airfoil  flow  field 
analysis,  has  been  implemented  in  the  computer  code.  The  analysis  of  Ref.  21 
follows  the  integral  turbulence  energy  procedure  of  Refs.  22  -  24  by 
utilizing  a  turbulence  function  a^,  where 


1  ^  1/2 
7  % 


(23) 


and  a^  is  taken  as  a  function  of  a  turbulence  Reynolds's  number  of  the  form 

•l  ■  •«  *.  -  *>]  <“> 

a  -  .0115 
o 

f(R^)  -  100  R^£  1  (25) 

f(R^)  -  68.1  R^  +  614.3  R^  >■  40 


10 


and  a  cubic  curve  was  fitted  for  values  of  between  1  and  40.  As 
previously  discussed,  Refs.  22  -  24  utilized  an  integral  form  of  the 
turbulence  kinetic  energy  and,  therefore,  R^  was  defined  as  an  average 
value.  g 

‘‘t  ’  T  ''t  /  7  ^o**  ''  ‘‘^s 

In  the  present  effort  was  defined  as  the  local  ratio  of  turbulent  to 
laminar  viscosity,  a^  was  evaluated  via  Eq.  24  and  Cy  related  to  a^  via 
Eq.  23. 


Two-Equation  Model  -  (k-c) 

Although  the  one-equation  approach  does  relieve  some  restrictions 
present  in  the  mixing  length  approach,  it  still  requires  specification  of  a 
length  scale.  The  k-e,  two-equation  turbulence  model  [20,  40]  in  which  both 
the  turbulence  kinetic  energy  and  the  turbulence  dissipation  rate  are 
governed  by  transport  equations  represents  a  more  general  model.  In  this 
approach  the  k-equation  is  as  given  in  Eq.  20,  but  the  algebraic  relation  for 
e  given  by  Eq.  22  is  replaced  by  the  following  transport  equation. 

(pe)  +  V»(puc)  ■  ^*('5^  D  :  ID  )  +  2uw.j(V^U)^-  c^P  (27) 


where  following  Reference  20 


and 


^  -  1.55 

[1  -  0.3  exp(-Rj,)^] 

2.0 


'2«» 


and  Rj  is  defined  as: 


»  - 

ejj 


However,  attempts  to  solve  Eqs.  20  and  27  without  modification  fail 
because  an  appropriate  boundary  condition  for  e  at  a  solid  boundary  is 
difficult  to  prescribe  such  that  Eq.  27  is  satisfied.  In  order  to  circumvent 
this  problem,  Eq.  27,  the  turbulence  dissipation  equation,  has  been  modified 
by  the  inclusion  of  the  term: 

-  2mu^  (v2u)2 
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which  following  Jones  and  Launder  [20]  is  included  in  order  to  obtain  better 
agreement  with  experimentally  determined  k  distributions  in  the  near-wall 
region.  In  addition,  Jones  and  Launder's  [20]  inclusion  of  the  term: 

-  2pv 

in  Eq.  28  is  a  numerical  device,  rather  than  a  physical  model,  which  allowed 
c>0  to  be  prescribed  as  a  function  boundary  condition  [20]. 

3(pk)  +  7  •  (pUk)  -  7  •  ^  7k  +  2W.J.  (  D  :  D  )  -  pe  -  2pv  (28) 
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COORDINATE  SYSTEM  DETAILS 


Equations  1-3  and  27  -  28  are  written  in  vector  form  in  which  the 
spatial  dependent  variables  represent  a  general  curvilinear  system  yJ 
(j  «  1,  2,  3).  In  principle  an  arbitrary  coordinate  system  can  be  used, 
however  experience  has  shotm  that  to  obtain  accurate  and  economical  numerical 
solutions  the  coordinate  system  used  must  meet  certain  criteria.  These 
include  the  ability  to  permit  accurate  implementation  of  boundary  conditions, 
sufficient  smoothness  of  coordinate  data  and  flexible  distribution  of 
computational  grid  points.  The  first  factor  to  be  considered  concerns  the 
treatment  of  boundary  conditions.  The  specification  of  boundary  conditions 
which  do  not  fall  upon  coordinate  lines  or  at  specific  grid  points  can 
present  a  difficult  problem  for  numerical  analyses.  Although  finite 
difference  molecules  can  be  constructed  for  boundary  points  which  do  not  fall 
upon  coordinate  lines,  such  molecules  may  have  an  unacceptably  high  spatial 
truncation  error  associated  with  them.  It  should  be  noted  that  while  the 
boundary  condition  problem  arises  in  both  viscous  and  inviscid  analyses,  it 
is  considerably  more  severe  in  viscous  flows  where  no-slip  conditions  on 
solid  walls  can  combine  with  boundary  condition  truncation  error  to  produce 
numerical  solutions  which  are  both  qualitatively  and  quantitatively  in 
error.  Thus,  the  first  property  required  of  any  coordinate  system  to  be  used 
in  a  viscous  analysis  should  be  that  boundary  surfaces  coincide  with 
coordinate  surfaces.  A  second  coordinate  system  requirement  is  sufficient 
smoothness  of  geometric  data.  This  requirement  is  more  stringent  than  simply 
requiring  the  existence  of  a  given  number  of  analytically  determined 
coordinate  transformation  derivatives;  the  coordinates  must  not  have  large 
changes  in  the  coordinate  geometry  metric  data  between  grid  points.  Such 
large  changes  can  cause  serious  deterioration  of  solution  accuracy  or  even 
prevent  obtaining  a  converged  solution.  The  final  item  to  be  considered  is 
the  need  to  obtain  high  grid  resolution  in  specified  regions  of  the  flow 
field.  For  example,  the  wall  region  of  viscous  flows  is  characterized  by 
rapid  changes  in  flow  variables,  therefore  the  computational  grid  should  have 
a  distribution  of  points  designed  to  resolve  these  gradients.  In  other 
regions  the  flow  variables  change  slowly  and  these  regions  should  have  a 
relatively  sparse  computational  grid  associated  with  them,  in  order  to 
maintain  an  economic  solution. 


A  recent  review  of  grid  generation  techniques  has  been  made  by  Thompson 
[25]  who  presents  an  extensive  overview  of  available  grid  generation 
techniques.  In  summary,  these  include  the  elliptic  methods  popularized  by 
Thompson  et  al  [26],  the  conformal  techniques  typified  by  Moretti's  work 
[27],  and  the  constructive  approach  such  as  that  of  Eiseman  [28]. 

The  governing  equations  for  the  present  problems  have  been  given  in  the 
previous  section  in  vector  form  (Eqs .  1-3  and  27  -  28).  However, 
implementation  of  a  solution  procedure  requires  that  these  equations  be 
transformed  into  an  appropriate  coordinate  system.  Therefore,  the  governing 
equations  written  in  a  cylindrical-polar  coordinate  system  are  transformed 
with  a  general  Jacobian  transformation  of  the  form 


y^  -  y^(xp  X2,  Xj,  t) 
T  «  t 


(29) 


where  (xj,  X2,  X3)  “  (r,  0,  z)  are  the  original  cylindrical  polar 
coordinates.  The  velocity  components  remain  the  components  (Up  U2, 

U3)  in  the  (xp  X2,  X3)  coordinate  directions,  respectively.  The  new 
independent  variables  yJ  are  the  computational  coordinates  in  the 
transformed  system.  The  coordinate  system  requirements  for  the  current 
application  may  be  represented  by  a  subset  of  the  general  transformation, 
Eq.  29, 


y^  »  y^(xp  Xj,  t)  (a) 

y^  -  y^(x2^ 

3  3  -  - 

y  -  y  (xj,  Xy  t)  (c) 

which  is  a  general  axisymmetric  time-dependent-trensformat ion.  For 
axisymmetric  flows,  Eq.  30b  reduces  to  y^  ■  X2»  and  all  derivatives 
3/3y2  are  assumed  to  be  zero.  The  transformation,  Eq.  30,  with  the 
axisymmetric  flow  assumption  has  been  utilized  in  the  computer  code  developed 
under  the  present  effort. 

Application  of  the  Jacobian  transformation  requires  expansion  of  the 
temporal  and  spacial  derivatives  using  the  chain  rule,  i.e.. 


(31) 
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and 


li  ■  ^  j 

3xj^  j*l  3y^ 


(32) 


where 


j  -  3 
y»i  =  -_ 
^  3x 


(33) 


The  relations  Eqs.  32  -  33  are  first  substituted  into  the  governing  equations 
written  in  cylindrical  polar  coordinates.  Then  the  resulting  equations  are 
multiplied  by  the  Jacobian  determinant  of  the  inverse  transformation, 


!5l 

3y^ 

3y^ 

3y3 

T  1 

3(xj,  Xj,  X3) 

(34) 

a/  1  2  3. 

3(y  ,  y  ,  y  ) 

3y^ 

3y^ 

3y^ 

!!! 

3y‘ 

3y^ 

3y^ 

and  the  equations  are  cast  into 

a  conservative  form  using 

the  following 

relations 

3 

3y\ 

Z 

-  0 

(35) 

j-l 

3y  j 

3j 

3 

3? 

t 

(36) 
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♦  1 
j-l 

-  -  0 

trhere 

^  »i 

>> 

III 

»i 

(37) 

;j 

^ 't 

=  Jyj 

•t 

It  should 

be 

noted  that  Equation  36  expresses  a  geometric 

conservation  law 

and  plays 

an 

important  role  in  enabling  the 

equations  of  motion  to  be  cast 

into  conservative  form. 
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The  particular  conservation  form  developed  implies  that  all  factors  involving 


/ 


:  I 

r 


the  radial  coordinate  r  *  remain  as  they  were  before  the  Jacobian 
transformation.  The  resulting  equations  are  presented  below  as  Eqs.  40  -  49. 
The  geometric  relations  Eq.  35  -  36  may  be  obtained  from  the 


••  i 

relations  for  y,^ 

and  ; 

“  J 

derivatives 

(e.g.. 

Reference 

-1 

y  ’1 

m 

"2.2 

*3.3  “ 

*2,3' 

*3,2 

-1 

_ 

y  *2 

m 

*3.2 

*1.3  ■ 

*3,3 

*1.2 

-1 

_ 

y  .3 

m 

*1.2 

*2.3  “ 

*1.3 

*2,2 

-2 
y  • 

m 

*2.3 

*3.1  ■ 

*2.1 

*3,3 

-2  ^ 

y  .2 

m 

*3.3 

*1.1  ■ 

*3.1 

*1.3 

-2 

_ 

y  .3 

m 

*1.3 

*2.1  ■ 

*1.1 

*2,3 

-3 

_ 

_ 

y  ‘1 

m 

*2.1 

*3.2  " 

*2.2 

*3,1 

-3 

y  .2 

m 

*3.1 

*1.2  ■ 

*3,2 

*1.1 

-3 

_ 

y  .3 

m 

*1.1 

*2.2  - 

*1,2 

*2,1 

?.t  - 

3 

P . 

-1 

k  - 

1 

(38) 


(39) 


The  transformed  governing  equations  may  be  written  in  the  following 
compact  form: 


3  •  .  3 


3(JW)  -  -  I  3  (Jy3,^W)  -  I  (8.  3  (JyJ,.F.) 

j  “  ^  3yJ  j  "  ^  3yJ 

+  Y.  _3_  (Jy^j  fj)  *  (V.i  2.)) 

3yj  3yJ 

+  JS  ♦  JC 
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(40) 


i. 


where 


3x. 

1 

Further,  the  coefficients  Y£,  Ci  are  given  by 


e-i.  0-1 

r  2  r  3 

Y  -  1 ,  Y  “  Y  •  1 

1  2  r  3 

?  ?  "i.  C  “1 

m  2  r  3 

r 


(41) 


(42) 


and  m  "  1  for  all  equations  except  the  X2  -  direction  momentum  equation, 
for  which  m  *  2.  The  vector  variables  used  in  Eq.  40  are  defined  as 


where  n  ■  1  for  i  ■  1  and  n  ■  0  for  i  ■  2,  3. 
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(46) 


(47) 


The  derivatives  required  in  Eqs.  46  -  47  must  be  expressed  in  terms  of  the 
computational  coordinates  yJ  using  the  chain  rule,  (Eq.  32). 

Finally,  the  vector  S  contains  source  terms  and  certain  differential 
terms  which  do  not  conform  to  the  basic  structure  of  Eq.  40,  and  the  vector  £ 
contains  the  additional  curvature  terms  due  to  the  cylindrical-polar 
coordinate  system. 


♦ 

S 


0 
0 
0 
0 

Dp  pe 

Dt 

U_  [20..  5.. 

T  ‘  ij  ij 

C,  e  [l'-(25.,  D.,)  ♦  2  WU..  (V^)‘ 
*  ^  * 


J  -pc  -  2pv  (7k^/^)^ 


P  £1  1 
^  k 


(48) 
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since  the  cases  described  in  this  report  have  widely  differing 
geometry,  specific  details  of  the  grid  transformations  will  be  postponed 
until  each  case  and  its  results  are  described  in  detail  in  subsequent 
sections . 

Solution  Procedure 

The  solution  procedure  used  to  solve  the  governing  equations  detailed 
previously  is  the  linearized  block  implicit  (LSI)  method  of  Briley  and 
McDonald  [15,  16],  which  has  been  applied  to  numerous  other  problems.  A 
detailed  explanation  of  this  method  can  be  found  in  the  work  of  Briley  and 
McDonald  [15,  16],  Repetition  of  this  is  unwarranted  here;  instead  a  summary 
of  the  method  can  be  found  in  the  appendix  to  this  report.  This  appendix 
also  contains  sections  dealing  with  linearization  and  time  differencing, 
treatment  of  diffusive  terms  and  the  splitting  scheme  of  the  LBI  method. 
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DISCUSSION  OF  RESULTS 


Under  the  present  effort  a  variety  of  test  cases  have  been  performed  in 
order  to  investigate  various  aspects  of  both  steady  and  unsteady  shock  wave¬ 
boundary  layer  interactions.  These  cases,  which  are  sumnarized  below,  are 
discussed  in  detail  in  this  section.  For  clarity,  geometry  descriptions  and 
boundary  condition  specification  are  described  subsequently  for  each  case. 

The  cases  presented  are  axisymmetric  unless  otherwise  stated,  and  are  as 
follows : 

(i)  A  normal  shock  wave-boundary  layer  interaction  in  a  tube 
(M»  *  1.44); 

(ii)  Steady  transonic  flow  over  a  bump  (Ma,  -  0.875); 

(iii)  Unsteady  transonic  flow  over  a  bump  (Mas  ==  0.875,  reduced 
frequency  =  0.175  and  0.004); 

(iv)  Flow  over  a  compression  ramp  (Mor>  -  2.0); 

(v)  Shock  impingement  on  a  flat  plate  boundary  layer  flow 
(M»  =  2.0). 

Normal  Shock  Wave-Boundary  Layer  Interaction  in  a  Constant  Area  Tube 

Normal  shock-wave-boundary  layer  interaction  is  a  frequently  occurring 
phenomenon  in  aeronautical  applications.  The  flow  field  contains  a  complex 
set  of  phenomena  including  shock  waves,  boundary  layers  in  adverse  pressure 
gradients  and  possibly  a  zone  of  separated  flow.  Even  in  cases  where  the 
shock  pattern  is  essentially  steady,  the  sudden  wall  static  pressure  rise 
accompanying  the  interaction  can  lead  to  significant  unbalanced  forces,  and 
the  high  heat  transfer  rates  typically  present  in  non-isothermal  flows  at 
flow  reattachment  can  lead  to  structural  problems.  The  first  case  considered 
simulates  a  normal  shock  wave-boundary  layer  interaction  occurring  at  modest 
supersonic  Hach  number  in  a  tube  of  constant  circular  cross  section.  The 
case  has  been  specified  to  coincide  with  the  experimental  data  of  Mateer, 
Brosh  and  Viegas  [30].  These  experimental  results  are  available  in 
sufficient  detail  to  enable  detailed  comparison  between  experiment  and 
prediction  to  be  performed. 

The  calculation  was  performed  for  an  inlet  Mach  number  of  1.44  with  an 
imposed  exit  to  inlet  pressure  ratio  of  2.0  and  a  Reynolds'  number  based  on 
upstream  boundary  layer  thickness  of  5.83  x  10^.  The  inlet  boundary  layer 
thickness  lAich  was  2.5  cm  was  also  used  as  the  reference  length  for  the 
calculation. 


Boundary  Conditions 

A  major  factor  in  obtaining  efficient  solutions  for  the  Navier-Stokes 
equations  is  specification  of  appropriate  boundary  conditions.  Boundary 
condition  specification,  which  has  been  the  subject  of  considerable  effort  in 
recent  years,  essentially  imposes  the  effect  of  the  environment  in  regions 
outside  the  computational  domain  on  the  problem  being  investigated. 
Specification  of  boundary  conditions  at  walls  is  relatively  straight  forward, 
but  proper  specification  at  inflow  and  outflow  boundaries  presents  a  more 
difficult  problem.  The  present  effort  follows  the  work  of  Briley  and 
McDonald  [33]  who  suggest  examining  the  characteristics  of  the  inviscid 
problem  for  guidance  at  inflow  and  outflow  boundaries.  Since  at  a  supersonic 
inflow  boundary  all  dependent  variables  can  be  specified  and  since  at  a 
subsonic  outflow  boundary  only  one  dependent  variable  can  be  specified, 
boundary  conditions  were  set  as  follows: 

(i)  Upstream  (supersonic  inflow)  boundary  - 

Function  boundary  conditions  for  each  dependent  variable 
(u,  w,  p,  h,  k  and  e)  with  profiles  consistent  with 
experimentally  determined  values  [30]. 

(ii)  Downstream  (subsonic  outflow)  boundary  - 

Second  derivative  of  u,  w,  h,  k  and  c  set  to  zero. 

Static  pressure  specified. 

(iii)  Wall  boundary  - 

No-slip  (i.e.  function  values  u«w“k"  e«0). 

Normal  Momentum  equation. 

Adiabatic  Wall. 

(iv)  Symmetry  Plane  - 

Function  value  u  ■  0. 

First  Derivative  of  w,  p,  h,  k  and  c  set  to  zero. 

A  second  major  consideration  in  obtaining  accurate  numerical 
simulations  of  flow  situations  is  specification  of  a  proper  computational 
grid.  Obviously,  code  size,  computer  core  size,  computer  run  costs,  etc. 
dictate  that  the  number  of  computational  grid  points  used  in  performing  a 
calculation  be  minimized.  However,  in  most  flow  situations  not  all  regions 
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of  the  flow  require  the  same  grid  resolution.  In  regions  containing  boundary 
layers  or  shock  waves  dependent  variables  change  rapidly  with  physical 
distance  and  these  flow  regions  require  high  grid  resolution.  Other  flow 
regions  may  be  satisfactorily  resolved  with  a  considerably  sparser  grid  point 
density.  As  previously  described,  the  governing  Navier-Stokes  equations  are 
written  in  a  general  coordinate  system.  However,  before  performing  a 
specific  calculation  a  coordinate  system  transformation  suitable  for  the 
problem  at  hand  must  be  provided.  Many  possible  methods  are  available  for 
grid  generation,  but  in  this  example  i^ere  both  boundary  geometry  and  flow 
structure  conform  to  well  defined,  geometric  shapes,  construction  of  the 
coordinate  system  is  preferable. 

For  the  axisymmetric  normal  shock  problem  of  this  section,  a  stretched 
Cartesian  grid  has  been  used  in  the  symmetry  plane.  For  this  purpose,  a 
modification  of  the  method  of  Roberts'  [311,  was  employed.  This  approach 
uses  a  sinh  function  for  the  purpose  of  clustering  points  in  the  interior  of 
a  domain,  to  resolve  the  shock  gradients,  and  a  tanh  function  at  boundaries 
to  resolve  boundary  layer  gradients.  A  typical  example  of  the  grid,  at  a 
particular  time,  is  given  in  Figure  1.  It  should  be  noted  in  this  connection 
that  care  must  be  exercised  in  using  the  sinh  function,  since  excessive 
stretching  can  lead  to  the  introduction  of  large  spatial  truncation  errors. 

Solution  Adaptive  Coordinate  System 

Specification  of  high  grid  resolution  in  the  vicinity  of  a  no-slip  wall 
is  obviously  required  to  resolve  a  boundary  layer.  However,  specification  of 
a  high  resolution  region  to  resolve  a  shock  is  not  as  simple  a  problem.  In 
some  flow  situations,  both  the  shock  location  and  angle  can  be  well  estimated 
and  in  these  cases  an  accurate,  a  priori  judgment  can  be  made  in  specifying 
the  region  of  high  resolution.  However,  in  many  cases  the  situation  is  not 
clear  because  the  shock  location  and/or  its  shape  may  not  be  known  a  priori. 
Since  adequate  grid  resolution  in  the  shock  region  is  required  to  obtain 
sharp  solutions  without  spurious  spatial  oscillations,  an  alternative 
procedure  must  be  used.  In  the  present  case,  a  normal  shock  was  expected  in 
the  tube.  Therefore,  an  orthogonal  system  with  one  set  of  coordinates 
parallel  to  the  wall  and  one  set  normal  to  the  wall  would  be  appropriate  if 
proper  regions  of  high  resolution  could  be  obtained  to  resolve  the  shock. 


With  the  basic  form  of  the  coordinate  system  established,  a  viable 
shock  tracking  strategy  was  developed  using  the  wall  pressure  gradient  and 
pressure  rises  to  determine  the  shock  location.  In  particular,  a  search  for 
the  maximum  pressure  gradient  location  was  used  to  establish  a  definition  for 
the  shock  center.  In  the  early  stages  of  calculation  a  "noisy"  solution 
often  exists.  Any  spurious  oscillations  in  the  solution  would  prove 
detrimental  to  accurate  shock  center  determination  since  large  gradients  are 
possible  and  could  conceivably  exceed  the  gradient  at  the  actual  shock 
center.  In  order  to  avoid  this  possibility,  a  "filtering  scheme"  was  applied 
to  the  process  of  shock  location.  Figure  2  demonstrates  diagrammt ically  the 
relevant  features  of  the  search  and  filtering  processes.  In  addition  to  the 
filtering  applied  to  the  adaptive  grid  shock  location  scheme,  the  movement  of 
the  mesh  was  controlled.  This  was  achieved  by  allowing  only  limited  mesh 
movement  at  each  time  step.  This  avoids  the  introduction  of  excessive 
temporal  truncation  error  which  could  slow  convergence. 

Essentially,  the  procedure  identifies  turning  points  in  the  wall  static 
pressure  distribution,  examines  the  change  in  pressure  between  subsequent 
turning  points  and  carries  out  a  search  for  the  maximum  pressure  gradient  in 
the  interval  having  largest  pressure  rise.  Sections  having  negative  pressure 
gradient  are  ignored  for  this  particular  problem.  This  approach  assumes  that 
the  shock  is  centered  at  the  point  of  maximum  streamwise  pressure  gradient  in 
the  interval  having  the  largest  positive  pressure  increase.  This  effectively 
filters  out  all  but  the  most  significant  pressure  rise  which  is  assumed  to  be 
the  rise  due  to  the  shock.  In  more  general  cases  where  multiple  shocks  occur 
or  where  separation  and  reattachment  give  rise  to  a  plateau  in  the  pressure 
distribution  a  less  restrictive  filter  would  be  required  if  resolution  of 
each  significant  pressure  rise  is  desired.  A  significant  pressure  rise  is  of 
course  problem  dependent,  therefore,  general  guidelines  are  difficult  to 
prescribe . 

Once  the  shock  center  was  located,  a  new  grid  was  constructed  by 
centering  the  sinh  function  at  the  new  shock  location.  The  transverse  grid 
was  not  changed  in  the  normal  shock  calculation  even  though  local  boundary 
layer  thickening  occurs  in  the  interaction  zone.  However,  the  transverse 
grid  transformation  was  constructed  to  take  account  of  the  expected 
thickening  and  thereby  ensure  adequate  boundary  layer  resolution  throughout 
the  calculation  domain.  The  entire  process  of  establishing  the  shock 


location  and  constructing  a  new  grid  was  performed  explicitly,  after  each 
time  step,  and  as  a  consequence  the  grid  motion  lags  the  shock  motion.  The 
grid  motion  was  accounted  for  in  the  governing  equations  through  inclusion  of 
terms  containing  derivatives  of  computational  coordinates  with  time. 

Results 

The  present  case  was  run  with  two  different  turbulence  models.  These 
were  a  mixing  length  model,  and  a  k  -  e  model;  both  of  which  have  been  de¬ 
scribed  earlier.  Before  beginning  the  calculation,  an  initial,  consistent 
approximation  to  the  flow  field  was  constructed.  This  was  obtained  by  assum¬ 
ing  the  Rankine-Hugoniot  relations  for  both  static  pressure  and  velocity 
ratios  across  the  shock,  at  an  assumed  streamwise  location,  and  applying  a 
tanh  function  to  blend  the  upstream  and  downstream  values  for  a  given  trans¬ 
verse  location.  By  also  assuming  constant  total  enthalpy  and  a  specified  up¬ 
stream  velocity  profile,  a  consistent  initial  flow  field  was  obtained.  The 
assumed  shock  location  did  not  take  into  account  the  influence  of  the  bound¬ 
ary  layer  on  shock  location  and  as  expected,  the  shock  location  began  to  move 
upon  initiation  of  the  calculation.  In  order  to  maintain  adequate  resolution 
in  all  phases  of  the  calculation,  the  adaptive  mesh  strategy  described  ear¬ 
lier  was  invoked.  The  calculation  was  first  performed  on  a  computational 
grid  with  41  transverse  and  31  streamwise  points,  using  the  mixing  length 
turbulence  model  described  earlier.  This  calculation  converged  within  150 
time  steps  with  an  overall  reduction  in  the  maximum  residual  in  the  field  of 
two  orders  of  magnitude.  The  residual  was  monotonically  decreasing  with 
increase  in  time  and  the  solution  was  stationary  when  the  calculation  was 
terminated . 

The  results  obtained  are  shown  in  Figs.  3  and  4  along  with  those 
obtained  using  a  two-equation,  k  -  c  turbulence  model.  In  the  current  work 
an  intermediate  calculation  was  also  performed  using  a  one-equation,  k  -  £ 
model.  This  additional  calculation  was  performed  in  order  to  obtain  a  suit¬ 
able  initial  k  field  for  use  in  starting  the  k  -  e  model  prediction. 

However,  more  recent  work  indicates  that  the  additional  step  is  not  necessary 
and  that  the  k  -  e  calculation  can  be  performed  by  using  the  mean  flow 
obtained  using  mixing  length  assumptions  and  using  initial  approximations  for 
k  and  e  determined  from  algebraic  relations.  Results  obtained  using  both 
mixing  length  and  k  -  e  models  are  shown  in  Pigs.  3-4  which  also  show 
comparison  with  experimental  data,  for  the  wall  static  pressure  through  the 


interaction  zone  and  velocity  profiles  at  the  streamwise  locations  indicated 
in  Fig.  4.  Figure  5  shows  the  shear  stress  profiles  evaluated  from  the  k  -  e 
results.  Closer  examination  in  Fig.  3,  vdiich  shows  the  wall  static  pressure 
distribution  in  the  interaction  zone,  reveals  good  agreement  with 
experimental  data  using  both  turbulence  models.  Figure  4  shows  velocity 
profile  plots  at  four  streamwise  locations.  Once  again,  agreement  between 
prediction  and  experiment  is  good  with  the  mixing  length  model  producing 
slightly  better  results. 

In  performing  this  calculation  the  artificial  dissipation  parameter  a, 
designed  to  remove  pre-  and  post-shock  induced  spurious  oscillations  and 
described  in  the  appendix,  initially  had  a  value  of  O.S.  Such  a  relatively 
large  value  certainly  smears  out  spurious  solution  oscillations,  but  can  also 
modify  the  solution  unrealistically.  Therefore,  after  convergence  with 
0  ~  O.S  its  value  was  reduced  to  0.05,  which  is  a  relatively  small  value,  in 
order  to  ensure  an  accurate  solution  free  from  arbitrary  smearing.  Values  in 
the  range  0.02S  <  o  <  O.I  have  been  used  by  the  authors  with  good  results  for 
a  variety  of  flow  problems  [42],  and  generally  a  further  reduction  below  o  > 
.OSleaves  the  solution  unaltered  with  eventually  the  spurious  oscillations 
returning  as  o  0. 

It  can  be  seen  that  mean  flow  quantities  are  well  predicted  by  the 
method,  but  apparent  discrepancies  exist  in  the  prediction  of  Reynolds  shear 
stress  profiles.  Fig.  5.  This  is  particularly  true  at  z/5  «  4  which  is  very 
close  to  the  interaction  zone.  It  is  interesting  to  note  that  the 
predictions  of  Mateer,  et  al  [30]  show  a  similar  disparity.  It  should  also 
be  noted  that  the  turbulence  measurements  were  obtained  using  hot  wire 
anemometry,  and  as  indicated  by  Mateer,  et  al  [30]  it  is  difficult  to 
determine  the  sensitivity  of  hot  wire  probes  to  velocity  and  density  changes 
in  transonic  flow.  The  discrepancy  occurs  in  the  region  of  the  shock 
boundary  layer  interaction  where  density  and  velocity  gradients  are  most 
significant.  As  a  consequence,  the  experimental  turbulence  data  must  be 
viewed  with  caution  in  this  area.  In  addition,  the  apparent  rapid 
dissipation  occurring  within  46«  of  the  interaction  region  seems 
unrealistic  and  inconsistent  with  known  turbulence  processes  (Fig.  5).  The 
problem  of  ascertaining  the  cause  of  the  discrepancy  would  be  difficult  to 
resolve  without  the  use  of  non-intrusive  turbulence  measuring  techniques. 
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To  gain  additional  insight  into  the  noraal  shock  wave-boundary  layer 
interaction  process,  it  is  informative  to  examine  contour  plots  of  static 
pressure  and  Mach  number.  Figure  6  shows  static  pressure  contour  plots  of 
the  whole  solution  domain  in  conjunction  with  a  coordinate  system  plot 
showing  the  relative  locations  of  grid  clustering  and  shock  location.  It  can 
be  seen  that  a  sharp  (shock)  pressure  rise,  extending  across  only  eight 
streamwise  grid  points  is  predicted  in  the  free  stream  region,  and  that  the 
effect  of  its  interaction  with  the  boundary  layer  is  to  diffuse  the  shock  and 
form  a  characteristic  lamda  foot  near  the  wall.  Figure  7  shows  the  static 
pressure  in  the  interaction  region  in  more  detail.  Figure  8,  which  presents 
velocity  vectors  in  the  shock  region  indicates  the  displacement  of  the 
boundary  layer  which  results  from  the  interaction  process.  It  should  be 
noted  that  the  present  calculation  strategy  ignores  the  possible  shock 
discontinuity  and  instead  continues  to  assume  the  solution  can  be  expanded  in 
a  Taylor  series  in  this  region.  The  present  procedure  'shock  captures' 
rather  than  'shock  fits'.  In  the  particular  flow  studied,  the  moderate 
super-sonic  Mach  number  encountered,  the  presence  of  the  boundary  layer,  and 
the  confined  nature  of  the  flow  all  help  ensure  that  in  reality  a  sharp 
discontinuity  in  the  flow  is  not  encountered,  thus  justifying  the  'shock 
capturing'  approach. 

The  normal  shock  boundary  layer  interaction  prediction  described  above 
demonstrates,  by  comparison  against  experimental  data,  an  ability  to 
accurately  and  economically  predict  transonic  shock  wave-boundary  layer 
interaction.  The  problem,  though  conceptually  simple,  provided  a  means  of 
validating  the  correct  operation  of  the  code  and  the  various  changes 
incorporated  to  allow  for  the  treatment  of  shock  waves.  The  great  similarity 
between  the  results  obtained  using  the  simple  mixing  length  and  the  k  -  e 
models  is  not  surprising.  This  is  because  the  mixing  length  was  originally 
developed  for  attached  boundary  layer  flows  emd,  at  least  in  regions  outside 
the  interaction  zone,  the  problem  considered  contains  discernable  boundary 
layer  characteristics.  Nevertheless,  the  problem  provided  a  useful  test  of 
the  k  -  e  model  for  which  the  mean  flow  results  show  equally  good  agreement 
with  data.  It  may  be  concluded  from  the  results  of  this  case  that  for 
related  normal  shock  wave-boundary  layer  interaction  problems  in  which  no 
flow  separation  occurs  use  of  a  mixing  length  model  will  provide  acceptable 
results. 
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Steady  Transonic  Flow  Over  an  Axisymmetric  Bum] 


The  second  portion  of  the  present  effort  focused  upon  the  effect  of 
longitudinal  curvature  on  the  shock  wave-boundary  layer  interaction  process. 
This  was  done  by  considering  the  problem  of  steady  transonic  flow  over  an 
axisyninetric  bump.  The  first  case  in  this  portion  of  the  study  was  based  on 
the  experiment  of  Johnson  and  Horstman  [32].  In  this  same  reference,  Johnson 
and  Horstman  also  report  the  details  of  their  own  calculation  for  the  same 
problem.  The  current  calculated  results  for  the  Johnson-Horstman  case  are 
compared  with  both  experiment  and  the  calculation  reported  in  Ref.  32.  The 
remaining  steady  transonic  bump  calculations  reported  are  variations  of  this 
original  case,  and  demonstrate  the  effect  of  increased  approach  boundary 
layer  thickness  and  of  longitudinal  wall  curvature  on  shock  inducement  and 
the  shock  wave-boundary  layer  interaction  process. 

Each  case  of  the  series  of  three  cases  was  performed  with  an  approach 
Mach  number  of  0.875  with  the  bump  mounted  on  a  cylinder.  The  experimental 
arrangement,  shown  in  Fig.  8,  used  a  hollow,  thin-walled  cylinder  to  avoid 
cylinder  leading  edge  effects  influencing  the  interaction  process.  There¬ 
fore,  in  simulating  this  experiment,  a  negligibly  thin  approach  boundary 
layer  was  specified  at  the  upstream  boundary.  The  bump  considered  is  a  cir¬ 
cular  arc  with  its  leading  edge  joined  smoothly  to  the  cylinder  surface  by  a 
circular  arc  fillit.  Using  the  chord  length  of  the  bump  as  the  reference 
length,  the  relative  dimensions  of  the  problem  are  indicated  in  Fig.  9.  The 
upstream  distance  from  the  bump  leading  edge  to  the  inflow  boundary  was  cho¬ 
sen  to  allow  development  of  a  boundary  layer  velocity  profile  consistent  with 
that  of  the  experiment.  The  other  steady  cases  completed,  although  not  com¬ 
pared  to  experiment,  are  compared  to  the  first  case  in  order  to  demonstrate 
in  general  terms  the  effects  of  increased  inlet  boundary  layer  thickness  and 
of  increasing  the  wall  curvature.  The  boundary  conditions  for  all  three 
cases  were  the  same  and  are  described  in  the  following  section. 

Boundary  Conditions 

Following  Briley  and  McDonald  [33]  analysis  of  the  characteristics  for 
the  bump  problem  shows  that  for  the  cases  of  interest  ^ere  the  inflow  and 
outflow  are  subsonic  and  the  supersonic  patch  expected  in  the  vicinity  of  the 
bump  does  not  extend  to  the  top  boundary  of  the  solution  domain  the  following 
boundary  conditions  are  appropriate. 
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(i)  Upstream  (subsonic  inflow)  boundary  - 

Function  condition  on  u-velocity  (transverse  velocity). 

Two  layer  model  on  w-velocity  and  static  enthalpy  [33]. 

See  following  text  for  a  brief  explanation. 

Second  derivative  of  static  pressure  set  to  zero. 

(ii)  Downstream  (subsonic  outflow)  boundary  - 

Second  derivative  of  u,  w  and  h  set  to  zero. 

Static  pressure  specified. 

(iii)  Wall  boundary  - 

No-slip  on  u  and  w. 

Normal  momentum  equation. 

Wall  temperature  specified. 

(iv)  Freestream  Boundary  - 

u  set  to  zero. 

First  derivative  of  w,  p  and  h  set  to  zero,  simulating  a  free 
slip-wall  condition. 

The  so-called  two-layer  model  [33]  used  at  the  inflow  boundary  is 
essentially  a  total  pressure  boundary  condition  applied  to  the  core  flow  with 
a  specified  boundary  layer  profile  shape  for  the  wall  region.  Matching  the 
static  pressure  at  the  edge,  defined  by  the  first  computational  point  from 
the  wall  at  which  |w|/|w|g,^  was  greater  than  or  equal  to  0.99  on  the 
previous  time  step,  enables  calculation  of  |w|  at  this  point.  This  provides 
the  required  normalizing  value  for  the  pre-specified  boundary  layer  profile 
shape.  Overall,  the  method  provides  a  mechanism  for  drawing  mass  flow  in 
order  to  satisfy  the  downstream  pressure  given  an  upstream  core  total  pres¬ 
sure.  This  specification  corresponds  to  the  usual  wind  tunnel  experiment 
where  stagnation  conditions  are  set  in  an  upstream  reservoir  and  static  pres¬ 
sure  is  set  at  some  downstream  location.  Specification  of  an  inlet  mass  flux 
could  be  accomplished  indirectly  by  varying  the  downstream  static  pressure, 
otherwise  attempts  to  specify  an  inflow  mass  flux  directly  leads  to  numerical 
problems . 

(i)  Large  Radius  of  Curvature  Bump;  Thin  Inlet  Boundary  Layer 
This  case,  for  t^ich  experimental  data  are  available,  simulates  flow  of 
a  developing  boundary  layer  over  a  moderately  thin  axisymmetric  bump.  Before 
beginning  the  calculation,  an  appropriate  initial  flow  field  was  con¬ 
structed.  This  was  achieved  by  imposing  a  linearly  varying  static  pressure 
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drop  betweea  inlet  and  exit  in  addition  to  an  assumed  upstream  velocity 
profile.  The  boundary  layer  thickness  upstream  was  set  to  a  value  of  0.0009c 
where  c  is  the  bump  chord  length  and  this  was  found  to  give  a  boundary  layer 
growth  matching  that  of  the  experiment  at  0.25  chords  ahead  of  the  bump 
leading  edge.  This  was  determined  by  comparing  the  computed  and  experi¬ 
mentally  reported  velocity  profiles  at  this  location.  The  details  of  the 
flow  field  were  then  obtained  by  assuming  flat  plate  flow  and  constant  total 
enthalpy  with  density  obtained  via  the  equation  of  state.  A  constructive 
technique  was  used  to  generate  a  computational  grid  of  31  transverse  and  51 
streamwise  points.  Later  the  number  of  grid  points  in  the  streamwise  direc¬ 
tion  was  reduced  to  A6  for  reasons  discussed  subsequently. 

The  distribution  of  streamwise  grid  points  for  this  problem  was  a 
compromise  between  the  desire  to  resolve  the  region  of  flow  acceleration  over 
the  bump  and  the  region  of  deceleration  through  the  shock  and  the  need  to 
maintain  an  economical  calculation.  In  view  of  this,  the  final  coordinate 
system  used  has  one  region  of  streamwise  grid  refinement.  This  region  spans 
the  bump  chord  and  has  a  minimum  grid  spacing  of  0.037c  at  the  trailing  edge, 
compared  with  a  local  maximum  grid  spacing  ahead  of  the  bump  of  0.523c,  and 
of  0.166c  aft  of  the  bump.  The  starting  flow  field  for  this  case  was  con¬ 
structed  by  assuming  the  flow  field  to  be  zero  pressure  gradient  flat  plate 
boundary  layer  flow.  In  other  words  a  self  similar  velocity  profile  was 
distributed  throughout  the  domain  with  density,  enthalpy  and  pressure 
specified  consistently.  No  attempt  was  made  to  introduce  an  initial  shock 
but  rather  the  flow  structure  was  allowed  to  develop  with  time.  Having 
established  a  suitable  starting  flow  field,  the  calculation  was  begun.  Solu¬ 
tion  convergence  was  achieved  within  about  150  time  steps  at  which  point  the 
maximum  residuals  of  the  problem  had  been  reduced  by  two  orders  of 
magnitude.  It  should  be  noted  that  the  calculation  at  this  time  was  changing 
about  some  mean  state  with  a  maximum  amplitude  of  about  IZ.  No  detailed 
study  of  accelerating  convergence  was  made  and  it  is  possible  that  con¬ 
vergence  could  be  obtained  with  significantly  fewer  time  steps.  In  the  final 
calculation  the  artificial  dissipation  parameter,  described  in  the  appendix, 
was  set  at  0.05.  The  calculation  was  performed  using  a  mixing  length  model 
to  account  for  turbulent  viscosity  variations  and  was  subject  to  the  boundary 
conditions  outlined  in  the  previous  section.  The  upstream  Mach  number  was 
0.875,  but  acceleration  over  the  curved  surface  of  the  bump  lead  to  locally 
supersonic  flow. 
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The  streamwise  wall  sCaCic  pressure  distribution  predicted  by  the 
current  method  is  presented  and  compared,  in  Fig.  10,  with  both  the 
experiment  and  calculations  of  Johnson,  Horstman  and  Bachalo  [32],  It  is 
imnediately  apparent  that  agreement  with  experiment  is  not  good  in  the 
interaction  region.  The  predicted  pressure  distribution  indicates  a  shock 
location  significantly  further  downstream  than  that  seen  experimentally. 
Similar  disparity  is  evident  in  the  results  of  Johnson,  et  al  vdio  attribute 
the  discrepancy  to  inadequacy  of  the  turbulence  models  used  in  their 
calculations.  Johnson,  et  al  [32]  used  two  different  models:  the  first 
being  the  Cebeci-Smith  [34]  mixing  length  model  and  the  second  the 
two-equation  model  of  Wilcox-Rubesin  [35].  In  view  of  the  close  agreement 
between  these  two  calculations  obtained  by  Johnson,  et  al  [32],  it  seems 
unlikely  that  turbulence  modeling  is  the  key  to  resolving  the  major 
discrepancy  between  calculation  and  experiment.  In  support  of  the  last 
statement,  and  purely  for  demonstration  purposes,  a  laminar  flow  calculation 
was  performed  using  the  method  of  this  report  and  subjected  to  the  same  inlet 
conditions,  boundary  conditions  and  pressure  ratio  of  the  turbulent  case. 

For  laminar  flow  at  the  same  Mach  number  as  the  turbulent  case  it  is  expected 
that  the  shock  would  adopt  a  position  closer  to  the  leading  edge  of  the 
bump.  Therefore,  a  laminar  calculation  should  indicate  the  possible 
sensitivity  of  the  calculation  to  the  turbulence  model.  Although  the  laminar 
calculation  indeed  moved  the  shock  forward  and  consequently  moved  the  results 
closer  to  experiment,  the  movement  was  not  sufficient  to  come  close  to 
matching  the  data.  Since  the  laminar  case  represents  an  extreme  perturbation 
in  the  viscosity,  and  since  the  two  calculations  of  Ref.  [32]  and  the 
turbulent  calculation  presented  here  are  essentially  in  agreement,  it  seems 
unlikely  that  physically  reasonable  perturbations  of  a  turbulence  model  would 
produce  better  agreement. 

More  likely  as  an  explanation  is  that  insufficient  experimental 
downstream  flow  and  geometry  details  were  provided,  leading  to  the 
application  of  incorrect  downstream  boundary  conditions  and  geometry  modeling 
in  both  the  Johnson-Horstman  calculation  and  that  reported  here.  It  should 
also  be  noted  that  Johnson  and  Horstman  [32]  applied  downstream  boundary 
conditions  at  a  distance  of  4.4  chords  from  the  leading  edge  in  order  to  be 
able  to  apply  first  derivative  boundary  conditions.  In  order  to  justifiably 
assume  that  flow  variables  were  not  changing  at  the  downstream  boundary,  they 
recalculated  the  flow  with  the  boundary  at  different  locations  and  imposed 
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Che  same  boundary  conditions.  However,  even  though  they  observed 
substantially  unchanged  results,  it  should  be  noted  that  Che  downstream 
boundary  of  their  calculation  was  placed  about  2  chords  further  downstream 
Chan  Che  diffusing  section  present  in  Che  experiment  and  as  a  result  did  not 
Cake  into  account  the  change  in  pressure  expected  in  the  transition  to  the 
diffuser,  shown  in  Fig.  9. 

In  Che  case  of  Che  current  calculation,  Che  downstream  boundary 
conditions  were  initially  applied  at  about  A  chords  downstream  of  the  leading 
edge  of  Che  bump  and  Che  magnitude  of  the  static  pressure  was  obtained  by 
extrapolation  from  Che  experimentally  determined  value  at  1.5  chords 
assuming,  as  Johnson  eC  al  implied,  Chat  Che  pressure  profile  was  flat  beyond 
Chat  location.  This  was  necessary  due  to  the  lack  of  experimental  data  at 
Che  end  of  Che  cylindrical  section,  approximately  2.5  chords  aft  of  Che 
leading  edge.  In  order  to  demonstrate  Che  importance  of  appropriate 
downstream  boundary  conditions,  two  additional  calculations  were  performed  in 
which  the  downstream  boundary  was  located  at  2.67  chords  aft  of  Che  leading 
edge.  In  this  case,  Che  number  of  sCreamwise  grid  points  was  cut  to  46. 

These  were  performed  with  the  same  boundary  conditions  as  described 
previously  but  with  Che  imposition  of  upstream  to  downstream  static  pressure 
ratios  of  1.0  and  1.05,  respectively. 

The  wall  static  pressure  for  both  of  these  cases  is  shown  superimposed 
on  Fig.  10.  It  can  be  seen  Chat  with  Che  downstream  boundary  located  at  Che 
start  of  Che  diffusing  section  of  the  experiment,  and  an  exit  static  pressure 
of  1.05,  Che  wall  pressure  profile  is  modified  and  Che  start  of  the  pressure 
rise  is  significantly  further  forward  Chan  in  the  case  with  Che  downstream 
boundary  at  4  chords  aft  of  Che  leading  edge.  With  a  pressure  ratio  of  1.05 
better  agreement  with  Che  experimental  data  is  achieved.  Even  though  Che 
minimum  and  maximum  pressures  are  greater  Chan  given  by  the  data,  the 
relative  shock  location  is  correct.  These  calculations  clearly  demonstrate 
Che  sensitivity  of  this  transonic  flow  configuration  Co  the  downstream 
condiC ions . 

Figures  11  and  12  show  constant  static  pressure  and  Mach  number 
contours  respectively  and  are  presented  in  order  to  demonstrate  Che  flow 
behavior.  Examination  of  Che  static  pressure  contours  in  Fig.  11,  shows  Che 
expected  expansion  as  Che  flow  accelerates  over  Che  bump,  followed 
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by  a  shock  as  the  flow  turns  back  parallel  to  the  cylinder  wall.  For  this 
case  flow  separation  is  of  limited  extent  with  the  separation  point  being  at 
z/c“0.69  while  reattachment  occurs  at  2/c“0.9.  Velocity  vectors  are  shown  in 
Fig.  13. 

An  additional  important  point  to  note  in  this  series  of  cases  is  the 
position  of  the  separation  and  reattachment  points  since  the  effective  wall 
shape  downstream  of  the  acceleration  determines  the  shock  location  which  in 
turn  affects  flow  separation  due  to  the  adverse  pressure  gradients  associated 
with  it. 

( ii)  Large  Radius  of  Curvature  Bump;  Thick  Inlet  Boundary  Layer 

This  second  axisymmetric  bump  case  was  performed  in  order  to 

demonstrate,  in  general  terms,  the  effect  of  increased  approach  boundary 
layer  thickness  on  the  transonic  shock  wave  boundary  layer  interaction.  The 
bump  geometry  used  was  identical  to  that  described  previously  and  the 
calculation  was  performed  with  the  imposition  of  an  inlet  to  exit  static 
pressure  ratio  of  1.0.  In  addition,  the  boundary  layer  thickness  imposed  at 
the  upstream  boundary  was  increased  from  0.0009c  to  0.049c.  Figures  14  and 
15  show  the  constant  static  pressure  and  Mach  number  contours  for  the  thick 
boundary  layer  case.  Comparing  Figs.  11  and  14  it  can  be  seen  that  the  shock 
in  Fig.  14  is  located  further  forward  than  that  in  Fig.  11,  and  that  it 
extends  further  into  the  free  stream.  The  thick  boundary  layer  shock  is 
weaker  and  the  shock  shape  closer  to  vertical,  thus  indicating  boundary  layer 
thickness  has  an  effect  on  strength  and  shape.  In  addition,  examination  of 
the  Mach  number  contours  reveals  an  extensive  separation  in  Fig.  15  which 
significantly  modifies  the  displacement  surface.  The  increased  boundary 
layer  thickness  associated  with  this  is  evident  in  Fig.  16,  which  shows 
velocity  vectors  at  selected  stations.  Comparisons  of  wall  static  pressure 
distributions  will  be  discussed  subsequently. 

(iii)  Small  Radius  of  Curvature  Bump;  Thin  Inlet  Boundary  Layer 

In  order  to  demonstrate  the  effect  of  bump  longitudinal  curvature  on 
shock  wave-boundary  layer  interaction,  one  final  steady  flow  case  was 
performed  in  which  the  thin  inlet  boundary  layer  of  the  original  case  was 
imposed  at  the  upstream  boundary  and  an  inlet  to  exit  static  pressure  ratio 
of  1.0  was  imposed.  These  conditions  were  imposed  on  a  cylinder/bump 


configuration  for  a  bump  having  a  radius  of  curvature  of  0.59c  compared  with 
1.35c  for  the  large  curvature  cases  above.  Figures  17  and  18  present 
constant  static  pressure  and  Mach  number  contours  respectively  while  Fig.  19 
shows  velocity  vectors  at  selected  streamwise  locations.  It  can  be  seen,  as 
expected  for  the  larger  bump,  that  a  pronounced  separation  occurs  at  the  bump 
trailing  edge.  In  addition,  while  a  shock  exists  it  does  not  extend  very  far 
into  the  core  flow.  However,  by  comparison  with  the  original  case,  having 
large  radius  of  curvature  and  thin  approach  boundary  layer,  the  start  of  the 
pressure  rise  is  located  significantly  closer  to  the  bump  leading  edge  at 
about  z  ~  0.7c.  Again  the  details  of  the  shock  pattern  including  strength 
and  shape  appear  sensitive  to  geometry. 

Steady  Transonic  Flow  Over  an  Axisymmetric  Bump:  Summary 

Among  the  main  observations  to  be  derived  from  the  foregoing  discussion 
of  results  are  the  following.  In  performing  calculations  of  this  type  it  is 
very  important  to  accurately  model  the  downstream  boundary  conditions. 

This  has  been  demonstrated  in  the  calculations  of  the  Johnson-Horstman  bump 
problem.  As  has  been  shown  in  Fig.  10,  specification  of  downstream  boundary 
conditions  appears  to  have  a  much  more  significant  impact  upon  results  than 
reasonable  (or  even  drastic)  changes  in  turbulence  model.  Also,  it  appears 
that  the  position  of  the  induced  shock  is  largely  determined  by  the  location 
and  extent  of  the  flow  separation,  which  in  turn  depends  upon  the  shock 
strength.  Thus,  an  analysis  containing  these  mutually  dependent  effects  is 
required  to  obtain  meaningful  predictions.  In  general,  both  increaseci 
approach  boundary  layer  thickness  and  more  severe  wall  curvature  tend  to  move 
the  separation  point  forward  since  even  without  the  shock  both  give  rise  to 
earlier  separation.  Increased  approach  boundary  layer  thickness  gives  rise 
to  premature  separation  because  of  the  boundary  layers  inability  to  sustain 
the  adverse  pressure  gradient  associated  with  the  flow  deceleration  through 
the  shock.  Increased  curvature  gives  rise  to  more  severe  acceleration  and 
deceleration  and  consequently  more  severe  adverse  pressure  gradients  which 
also  adversely  affects  the  ability  of  the  boundary  layer  to  avoid 
separation.  The  presence  of  the  shock  tAiich  appears  as  a  result  of  turning 
the  locally  supersonic  flow  enhances  the  pressure  gradients  giving  rise  to 
some  siodif icat ion  of  the  location  of  the  separation,  and  its  extent. 
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Figure  20  gives  e  comparison  of  the  wall  pressure  distribution  for  each 
of  the  above  cases.  It  is  apparent  chat  the  case  having  a  smaller  radius  of 
curvature  bump  yields  a  pressure  distribution  consistent  with  intensive  flow 
acceleration  over  Che  bump.  This  is  indicated  by  Che  significantly  lower- 
minimum  wall  pressure  achieved.  Both  this  case  and  Che  case  having  a  thick 
approach  boundary  layer  exhibit  flattening  of  the  pressure  rise  typical  of 
that  occuring  when  shock  induced  separation  occurs.  The  start  and  end  of 
this  "plateau"  coincides  in  these  calculations  with  Che  separation  and 
reattachment  points.  No  such  flattening  is  apparent  in  Che  larger  radius  of 
curvature,  Chin  boundary  layer  case  which  is  consistent  with  Che  lack  of 
significant  separation  in  this  case.  Considering  the  overall  wall  pressure 
distribution,  it  appears  Chat  Che  wall  pressure  is  more  sensitive  to  changes 
in  Che  bump  curvature  Chan  Co  changes  in  Che  boundary  layer  thickness.  This 
indicates  bump  curvature  may  be  an  extremely  important  parameter  in 
axisymmetric  interactions. 

Unsteady  Transonic  Flow  Over  an  Axisymmetric  Bump 

Unsteady  transonic  flow  is  a  frequent  occurence  in  aeronautical 
engineering  applications,  and  consequently  prediction  of  such  phenomena  are 
of  interest.  Examples  of  unsteady  transonic  flow  are  the  helicopter  rotor 
field,  Che  surge  phenomenon  occurring  in  transonic  turbomachinery  flows  and 
flight  control  surface  buzz.  Therefore,  under  Che  present  effort,  cases  of 
unsteady  transonic  flow  were  considered  with  Che  major  aims  being  Co 
demonstrate  unsteady  transonic  flow  calculations  and  consider  both  a 
quasi-steady  and  time-dependent  flow  configuration  using  the  time-dependent 
Navier  Stokes  equations.  The  specific  flow  configuration  considered  was 
obtained  by  imposing  a  periodically  perturbed  exit  static  pressure  on 
transonic  flow  over  an  axisymmetric  btmip.  The  bump  geometry  modeled  is 
identical  Co  that  described  earlier  in  which  a  large  radius  of  curvature  bump 
was  mounted  on  an  axisymmetric  cylinder,  with  negligibly  thin  approach 
boundary  layer. 

The  boundary  conditions  for  this  problem  were  described  earlier  with 
Che  exception  of  Che  exit  static  pressure  boundary  condition.  The  exit 
static  pressure  was  applied  as  a  time-dependent  function  of  the  form: 

p(t)  ■  p  &p  .  cos  (ut) 
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where  p  is  the  mean  exit  pressure,  Ap  is  the  perturbation  amplitude,  u)  is  the 
frequency  and  t  is  time.  Rather  than  use  (s,  it  is  convenient  to  introduce 
the  reduced  frequency  f^  based  upon  boundary  layer  scales  defined 
as  u)5ga/Uas  where  is  the  boundary  layer  thickness  at  the  leading  edge 
of  the  bump  and  UU  the  edge  velocity  at  this  same  location. 

For  the  two  cases  described  here,  p,  ~  1.05,  Ap  -  0.05  and 
fr  -  0.175  and  0.004. 

(i)  High  Frequency  Case  (fr  “  0.175) 

The  unsteady  flow  cases  described  here  were  performed  using  a  mixing 
length  model  to  account  for  variations  in  turbulent  viscosity.  Use  of  such  a 
simple  model  seems  reasonable  in  this*  case  since  the  main  aim  is  the  econom¬ 
ical  demonstration  of  the  capability  of  the  present  method  to  compute 
periodically  unsteady  flows  and  to  examine  the  effect  of  frequency  upon  the 
flow.  Furthermore,  appropriate  time-dependent  multi-equation  models  are  not 
yet  %fell  established.  The  starting  flow  field  for  this  case  was  the  con¬ 
verged  steady  case  described  previously  in  which  the  pressure  ratio  was  set 
to  1.0,  but  as  noted  in  the  introduction  to  this  section,  the  time  mean  pres¬ 
sure  for  these  cases  was  set  to  1.05.  However,  starting  the  calculation  such 
that  the  initial  exit  pressure  for  the  case  was  1.0  assured  a  gradual  transi¬ 
tion  from  steady  flow  to  unsteady  flow  calculation.  In  other  words  the 
calculation  was  begun  at  a  turning  point  in  the  exit  pressure  temporal  varia¬ 
tion;  i.e.,  dpg/dt"0.  Examination  of  the  time  history  of  the  wall  static 
pressure,  given  in  Fig.  21,  reveals  the  frequency  response  of  the  flow  to 
cyclic  variation  in  the  imposed  pressure  ratio.  Shown  in  the  figure  are  the 
static  pressure  variations  at  the  three  locations  indicated.  These  locations 
are  at  the  downstream  boundary,  at  z/c  *  0.75,  the  steady  flow  shock  location 
and  at  z/c  ■  0.61  the  experimentally  determined  steady  flow  shock  location. 
The  latter  two  locations  were  chosen  for  convenience  and  are  not,  of  them¬ 
selves,  of  major  physical  significance.  The  initial  wall  pressure  transient 
is  a  response  to  the  newly  imposed  mean  pressure  ratio  and  takes  7w  radians 
before  the  first  turning  point  is  reached.  Between  7Tr  radians  and  lOx 
radians  the  response  of  the  wall  pressure  at  z/c  ■  .75  developed  a  periodic 
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oscillation  corresponding  to  the  sinusoidal  oscillation  of  the  back 
pressure.  Examination  of  the  response  at  the  start  of  the  pressure  rise  at 
z/c  ■  .61  shows  no  appreciable  variation  once  the  initial  transient  adjusting 
the  mean  pressure  is  over.  The  implication  is  that  after  adjusting  to  the 
new  mean  value,  the  shock  location  is  essentially  unaffected  by  downstream 
oscillations  at  this  frequency.  That  is  at  this  relatively  high  frequency  of 
oscillation  disturbances  propagated  upstream  from  the  perturbed  exit  pressure 
are  attenuated  by  the  shock  and  as  such  do  not  affect  the  shock  location. 
Comparison  of  this  result,  at  the  instant  at  which  the  pressure  ratio  is 
about  1.05,  with  the  steady  state  solution  having  a  pressure  ratio  of  1.05 
shows  no  appreciable  difference  in  the  wall  static  pressure  profile  in  the 
vicinity  of  the  shock. 

(ii)  Low  Frequency  Case  (f^  “  0.004) 

In  order  to  show  the  flow  phenomena  for  a  much  reduced  exit  forcing 
frequency  a  low  frequency  case  was  initiated.  This  calculation  was  begun  by 
restarting  from  the  previous  high  frequency  case  but  with  the  new  frequency 
imposed  on  the  exit  pressure  variation.  For  this  case,  as  in  the  previous 
case,  the  mean  static  pressure  at  exit  was  1.05  and  the  amplitude  of  the 
oscillation  was  0.05.  However  the  reduced  frequency  of  oscillation  was 
reduced  to  0.004  in  order  to  attempt  to  appreciably  influence  the  location  of 
the  shock.  Modifying  the  frequency  of  oscillation  instantaneously  introduces 
a  transient  into  the  calculation,  which  was  allowed  to  develop  over  several 
cycles  before  a  periodic  response  to  the  oscillating  exit  pressure  was 
detected.  Two  additional  cycles  were  then  computed  in  order  to  verify 
periodicity.  Figure  22  shows  the  wall  static  pressure  variation  at  a  number 
of  streamwise  locations  including  the  exit  plane.  This  figure  indicates 
periodic  response  to  the  exit  forcing  frequency,  but  Fourier  analysis  of  the 
results  also  reveals  that  the  response  in  the  vicinity  of  the  shock  contains 
higher  harmonics  which  tend  to  distort  the  response  locally.  In  support  of 
this.  Figure  23  shows  the  streamwise  variation  in  anplitude  of  the  primary 
and  first  harmonics  of  the  wall  static  pressure  response.  It  can  be  seen 
that  far  upstream  and  downstream  the  response  remains  almost  purely 
sinusoidal  since  the  first,  and  higher  harmonics  are  negligible.  However  in 
the  vicinity  of  the  shock  the  amplitude  of  both  the  primary  and  first 
harmonics  peak  dramatically.  This  is  an  expected  result  of  the  motion  of  the 


37 


shock  since  in  the  vicinity  of  the  start  of  the  shock  the  time  mean  pressure 
profile  exhibits  a  minimum;  and  with  the  shock  oscillating  about  its  time 
mean  location,  wall  points  in  this  vicinity  (z/c  “  0.61)  experience  a  wider 
range  of  pressure  variation  than  those  upstream  and  downstream  of  the  bump 
where  the  time  mean  pressure  profile  is  a  relatively  slowly  varying  function 
of  z/c.  It  is  also  interesting  to  note  that  the  time  mean  pressure  profile 
is  indistinguishable  from  the  steady  flow  pressure  profile  having  an  exit 
pressure  of  1.05.  Figure  24  shows  a  sequence  of  static  pressure  contour 
plots  for  this  problem  and  indicates  the  movement  of  the  shock  as  a  function 
of  time. 


Unsteady  Flow  Cases;  Summary 

The  unsteady,  transonic  flow  calculations  performed  and  reported  in  the 
previous  section  are  for  two  reduced  frequencies  in  the  high  and  low  range 
respectively.  It  is  observed  that  for  high  frequency  of  exit  pressure  oscil¬ 
lation  modification  to  the  flow  field  are  confined  to  the  vicinity  of  the 
downstream  boundary.  The  shock  location  is  unaffected  by  the  perturbations. 
In  contrast  the  low  frequency  case,  having  the  same  forcing  function  anpli- 
tude,  exhibits  a  periodic  response  to  the  downstream  pressure  fluctuations. 
This  case  in  which  the  shock  moves  periodically  about  its  time  mean  location 
produces  large  (by  comparison  with  the  input)  amplitude  fluctuations  in  the 
vicinity  of  the  shock's  time  mean  location. 

Steady  Flow  Over  An  Axisymmetr ic  Compression  Ramp 

Many  %K>rkers  have  addressed  the  problem  of  calculating  oblique  shock 
wave  inducement  at  two-dimensional  compression  corners  in  supersonic  approach 
flow  [9].  However,  the  axisymmetric  compression  corner  problem  has  received 
little  attention.  Therefore,  as  part  of  the  present  effort  consideration  was 
given  to  the  problem  of  axisymnetric  compression  corners.  Before  performing 
the  calculation,  which  in  this  case  was  for  an  approach  flow  having  a  Mach 
number  of  2.0  and  a  23*  axisymmetric  compression  ramp,  a  suitable  starting 
flow  field  and  coodinate  system  were  required.  In  constructing  the 
coordinate  system  use  was  made  of  the  known  exact  incompressible  potential 
flow  solution  for  flow  over  an  axisymmetric  ramp  in  order  to  generate  the 
streamwise  coordinate  lines.  Adopting  the  streamlines  obtained  from  the 
potential  flow  solution  and  applying  an  appropriate  grid  stretch  in  the 
transverse  direction  to  resolve  the  boundary  layer  region  provided  one  family 


of  coordinate  lines.  The  transverse  coordinate  lines  were  obtained  by 
defining  the  end  points  of  an  assumed  oblique  shock  at  an  angleextracted  from 
the  inviscid  shock  wave  relations,  and  using  these  points  on  the  top  and 
bottom  boundaries  as  the  centering  location  about  which  the  remaining 
strearawise  grid  points  are  clustered.  The  second  family  of  coordinate  lines, 
therefore  consisted  of  a  series  of  straight  lines  joining  points  on  the  top 
and  bottom  boundaries  of  the  computational  domain.  Fig.  25.  These  lines  are 
vertical  at  the  upstream  boundary  but  as  a  result  of  applying  different  grid 
stretching  on  the  top  and  bottom  boundaries  the  angle  each  line  makes  with 
the  horizontal  varies  throughout  the  computational  domain  such  that,  in  the 
vicinity  of  the  expected  shock,  the  grid  lines  align  with  the  shock 
direction.  At  the  downstream  boundary  coordinate  lines  are  once  again 
vertical.  Having  generated  the  grid  the  initial  velocity  field  was  aligned 
in  the  direction  of  these  coordinate  lines.  To  account  for  the  expected 
shock  the  initial  flow  field  was  constructed  by  considering  the  flow  field  to 
consist  of  three  zones.  In  the  first,  upstream,  zone  pre-shock  velocities 
and  pressures  where  applied.  Downstream  post  shock  values  where  obtained 
from  the  shock  relations,  and  in  the  intermediate  region  upstream  and 
downstream  values  were  smoothly  blended  to  avoid  an  initial  discontinuous 
representation  of  the  shock.  Having  determined  both  the  velocity  and  static 
pressure  fields,  the  assumption  of  constant  total  enthalpy  was  invoked  to 
enable  the  density  and  temperature  fields  to  be  determined. 

The  resulting  flow  field  constituted  a  reasonable  approximation 
for  use  as  the  initial  flow  field  for  the  problem.  The  calculation 
initiated  from  this  flow  field  converged  within  150  time  steps;  although  in 
developing  the  strategy  for  running  the  case  approximately  500  time  steps 
were  performed,  it  is  estimated  that  to  rerun  the  same  or  a  similar 
calculation  with  all  the  appropriate  changes  incorporated  150  time  steps  is  a 
realistic  estimate  for  convergence.  At  convergence,  the  maximum  residuals 
for  the  problem  had  been  reduced  by  two  orders  of  magnitude,  at  which  point 
the  maximum  change  in  the  solution  on  subsequent  time  steps  was  less  than  1 
per  cent,  and  was  continually  diminishing  with  time. 

The  boundary  conditions  used  for  the  wedge  problem  were  based  upon  a 
characteristic  analysis  of  the  inviscid  flow  equations  which  indicates  that 
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for  Che  equations  used  four  physical  conditions  are  appropriate  at  Che 
upstream  boundary.  The  conditions  specified  for  the  problem  were: 

(i)  Upstream  (supersonic  inflow)  boundary  - 

Function  conditions  on  u-velocity  and  w-velocity  with  a 
bounday  layer  profile  specified. 

Function  conditions  on  density  and  enthalpy. 

(ii)  Downstream  (supersonic  outflow)  boundary  - 

All  dependent  variables  extrapolated  by  imposing  zero  second 
derivatives . 

(iii)  Wall  boundary  - 

No-slip  on  u-  and  w-velocities . 

Normal  derivatives  of  static  pressure  and  enthalpy  set  Co  zero, 

(iv)  Free  stream  boundary  - 

First  derivative  in  the  coordinate  direction  set  to  zero  for  all 
dependent  variables. 

The  assumption  made  in  applying  the  above  free  stream  boundary  conditon 
is  Chat  ahead  and  after  Che  shock,  Che  free  stream  is  expected  to  be  inviscid 
and  approximate  slug  flow.  In  addition,  in  Che  vicinity  of  the  shock  flow 
variables  in  the  direction  along  Che  shock  do  not  vary.  Therefore,  by 
arranging  Che  coordinate  system  Co  align  with  Che  shock  enables  extrapolation 
from  the  field  to  Che  free  stream  boundary  in  Che  coordinate  direction  to  be 
a  reasonable  boundary  condition. 

The  results  presented  are  for  a  calculation  in  which  the  mixing  length 
model  described  earlier  has  been  used  to  represent  Che  effect  of  turbulence 
on  Che  mean  flow  development  of  supersonic  flow  over  an  axisymmeCric 
compression  ramp.  These  results  are  sho«m  in  Figs.  26  to  28  and  include  wall 
and  freestream  static  pressure  distribution,  constant  density  and  constant 
static  pressure  contours.  Examination  of  Figs.  26  and  27  show  a  sharp  shock 
located  at  Che  corner.  The  upstream  influence  of  this  corner  is,  as 
expected,  very  limited  with  a  marginal  increase  in  pressure  ahead  of  Che 
corner . 

The  wall  static  pressure  distribution  of  Fig.  28  shows  a  very  sharp 
pressure  rise  Co  a  peak  pressure  of  2.65.  The  rise  in  pressure  does  not 
exhibit  Che  plateau  often  expected  in  induced  shock  flows.  In  general,  such 
a  plateau  is  indicative  of  a  separation  region.  In  this  case,  a  very  small 
separation  zone  exists  with  separation  and  reattachment  at  z  ■  -0.23  and 
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0.023,  respectively.  At  reattachment,  the  pressure  begins  to  fall  smoothly 
and  monotonically  down  to  a  level  of  2.35.  In  contrast  to  the  sharp  shock 
indicated  by  the  wall  static  pressure  distribution  the  pressure  distribution 
at  the  freestream  boundary  shows  a  steady  rise  in  pressure  through  a 
relatively  diffuse  shock  with  a  maximum  pressure  of  2.2.  This  pressure  rise 
corresponds  to  an  expected  inviscid  pressure  rise  of  1.2k.  These  results 
were  obtained  with  the  artificial  dissipation  parameter  o  set  to  0.1  which 
experience  has  shown  to  be  a  low  value.  Significantly,  the  results  are  free 
of  spurious  oscillation  in  spite  of  this  low  artificial  dissipation. 
Experience  with  other  similar  codes  has  shown  o  -  0.1  to  produce  acceptable 
solutions  tdtich  compare  well  with  data. 

In  the  case  described  here  the  lack  of  experimental  data  prevents 
comparisons  but  in  general  it  can  be  stated  that  the  solutions  appear 
physically  plausible. 

Oblique  Shock-Wave  Impingement  on  a  Flat  Plate  Boundary  Layer 

While  the  supersonic  flat  plate  boundary  layer  shock  impingement  problem 
has  been  addressed  by  a  number  of  workers  employing  Navier-Stokes  analyses, 
e.g.  [36],  most  have  used  a  Cartesian  coordinate  system  on  which  to  perform 
their  calculation.  The  main  criticism  of  these  calculations  centers  around 
the  shock  capturing  used  to  treat  the  impinging  and  reflected  shocks.  The 
use  of  a  Cartesian  grid  to  predict  an  oblique  shock  can  result  in  spatial 
truncation  error  problems  as  a  result  of  poorly  resolving  the  shock.  In 
order  to  alleviate  such  problems,  a  coordinate  system  was  devised  which 
aligns  with  both  the  impinging  and  reflected  shocks  thereby  maintaining 
adequate  resolution  of  both  shocks  with  a  limited  number  of  grid  points, 

Fig.  30. 

A  calculation  was  performed  with  this  coordinate  system  for  a  freestream 
Mach  number  M  "  2.96  and  an  impinging  shock  wave  created  by  a  10.89*  shock 
generator.  The  Reynolds  number  was  sufficiently  high  so  that  the  flow  was 
turbulent,  and  a  mixing  length  turbulence  model  was  employed.  The  coordinate 
system  shown  in  Fig.  30  was  allowed  to  adapt  to  the  impinging  and  reflected 
shocks  so  that  very  sharp  shock  definition  was  obtained,  as  shown  by  the 
pressure  contours  in  Fig.  31.  Even  so  the  predicted  pressure  variation  along 
the  wall  did  not  compare  very  well  with  the  available  experimental  data  for 
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this  case  [36).  In  particular,  the  predicted  wall  pressure  showed  a  much 
more  rapid  rise  than  that  observed  experimentally  with  no  plateau  which  is 
character  it ic  of  a  significant  separation  region.  While  the  calculation  did 
predict  separation,  its  length  and  extent  was  not  sufficiently  large  to 
modify  the  wall  static  pressure  distribution  significantly. 

In  the  calculation  an  artificial  dissipation  parameter  of  o  >  0.5  was 
employed,  and  all  attempts  to  reduce  this  value  were  unsuccessful.  It  was 
not  known  whether  the  coordinate  system,  the  artificial  dissipation  or  the 
relatively  simple  mixing  length  turbulence  model  was  responsible  for  the 
inaccuracies  in  the  predictions.  Therefore,  a  study  was  undertaken  to 
identify  the  cause  of  the  discrepancies.  First,  an  oblique  shock  impinging 
on  a  flat  plate  laminar  boundary  layer  was  considered  using  a  mesh  similar  to 
that  shown  in  Fig.  30.  The  case  chosen  had  a  freestream  Mach  number  M  =  2.0, 
and  a  shock  impinging  at  an  angle  of  about  32.6  degrees.  As  in  the  previous 
case,  the  predicted  wall  pressure  rise  was  more  rapid  than  that  observed 
induced  separation  region  was  not  as  large  as  that  indicated  by  the 
experimental  skin  friction  measurements.  Also,  in  this  case  it  was  not 
possible  to  reduce  the  artificial  dissipation  parameter  below  a  value  of 
o  »  0.5. 

Therefore,  a  calculation  was  performed  for  the  shock  impinging  on  a 
laminar  boundary  layer  using  a  Cartesian  grid  with  50  transverse  points  and 
100  streamwise  points.  The  computational  domain  extended  from  0.6  inch  to 
3.1  inches  downstream  of  the  plate  leading  edge,  and  from  0.0  to  2.0  inches 
normal  to  the  wall.  The  inflow  boundary  layer  thickness  (.02  inch)  for  a 
fully  developed  turbulent  profile  was  chosen  such  that  the  boundary  layer 
velocity  profile  was  reproduced  reasonably  well  at  the  first  data  measurement 
station,  1.1  inches  from  the  plate  leading  edge.  In  this  case  the  shock  wave 
impinges  at  about  1.96  inches,  and  the  shock  has  little  effect  on  the 
velocity  profile  at  the  first  measurement  station.  The  wall  pressure 
distribution  for  this  case  is  shown  in  Fig.  32,  where  it  is  seen  that  initial 
pressure  rise  begins  somewhat  upstream  of  the  experimental  rise  [41]  and  a 
plateau  longer  than  that  observed  experimentally  is  predicted.  The 
calculation  also  yielded  a  separation  region  extending  from  about  1.5  to  2.5 
inches  which  is  about  0.5  inches  longer  than  the  experimental  separation 
region  as  determined  from  skin  friction  meaurements.  This  larger  calculated 
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separation  region  is  consistent  with  the  longer  wall  pressure  plateau  shown 
in  Fig.  32. 

In  this  Cartesian  grid  calculation  an  artificial  dissipation  parameter 
of  o  >  0.5  was  required,  and  a  smooth  solution  could  not  be  obtained  with  a 
smaller  value.  As  a  result,  the  predicted  incident  and  reflected  shock  waves 
were  excessively  smeared  even  though  the  calculation  utilized  a  large  number 
of  mesh  points.  Thus,  the  pressure  rise  across  the  shock  extended  over  a 
larger  physical  region  beginning  further  upstream  than  experimentally 
observed,  but  with  a  predicted  pressure  gradient  large  enough  to  cause 
separation  of  the  laminar  boundary  layer. 

In  conclusion,  the  Cartesian  grid  calculation  of  wall  pressure  for  a 
shock  wave  impinging  on  a  laminar  boundary  layer  yielded  reasonable  results 
with  quantitative  discrepancies  attributable  to  excessive  shock  smearing  due 
to  artificial  dissipation.  Modifications  to  the  form  and  magnitude  of  the 
artificial  dissipation  required  for  shock  capturing  are  necessary  to  reduce 
these  errors.  Also,  shock  adapting  meshes  such  as  that  shown  in  Fig.  30 
should  be  further  investigated  because  of  the  potential  for  sharp  shock 
predictions  with  a  minimum  number  of  mesh  points. 
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APPENDIX  -  SOLUTION  PROCEDURE  [17] 
Background 


The  solution  procedure  employs  a  consistently-split  linearized  block 
implicit  (LBI)  algorithm  which  has  been  discussed  in  detail  in  [15,  16]. 

There  are  two  important  elements  of  this  method; 

(1)  the  use  of  a  noniterative  formal  time  linearization  to 
produce  a  fully-coupled  linear  multidimensional  scheme  which 
is  written  in  "block  implicit"  form;  and 

(2)  solution  of  this  linearized  coupled  scheme  using  a  consistent 
"splitting"  (ADI  scheme)  patterned  after  the  Douglas-Gunn  [37] 
treatment  of  scalar  ADI  schemes. 

The  method  is  thus  referred  to  as  a  split  linearized  block  implicit  (LBI) 
scheme.  The  method  has  several  attributes: 

(1)  the  noniterative  liearization  is  efficient; 

(2)  the  fully-coupled  linearized  algorithm  eliminates  instabilities 
and/or  extremely  slow  convergence  rates  often  attributed  to  methods 
which  employ  ad  hoc  decoupling  and  linearization  assumptions  to 
identify  nonlinear  coefficients  which  are  then  treated  by  lag  and 
update  techniques; 

(3)  the  splitting  or  ADI  technique  produces  an  efficient  algorithm 
which  is  stable  for  large  time  steps  and  also  provides  a  means  for 
convergence  acceleration  for  further  efficiency  in  computing  steady 
solut ions ; 

(4)  intermediate  steps  of  the  splitting  are  consistent  with  the 
governing  equations,  and  this  means  that  the  "physical"  boundary 
conditions  can  be  used  for  the  intermediate  solutions.  Other 
splittings  which  are  inconsistent  can  have  several  difficulties  in 
satisfying  physical  boundary  conditions  [16]. 

(5)  the  convergence  rate  and  overall  efficiency  of  the  algorithm  are 
much  less  sensitive  to  mesh  refinement  and  redistribution  than 
algorithms  based  on  explicit  schemes  or  which  employ  ad  hoc 
decoupling  and  linearization  assumptions.  This  is  important  for 
accuracy  and  for  computing  turbulent  flows  with  viscous  sublayer 
resolution;  and 
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(6)  the  method  is  general  and  is  specifically  designed  for  the 

complex  systems  of  equations  which  govern  multiscale  viscous  flow 
in  complicated  geometries. 

This  same  algorithm  was  later  considered  by  Beam  and  Warming  [38],  but  the 
ADI  splitting  was  derived  by  approximate  factorization  instead  of  the 
Douglas-Gunn  procedure.  They  refer  to  the  algorithm  as  a  "delta  form" 
approximate  factorization  scheme.  This  scheme  replaced  an  earlier  non-delta 
form  scheme  [39],  which  has  inconsistent  intermediate  steps. 

Split  LBI  Algorithm 

Linearization  and  Time  Differencing 

The  system  of  governing  equations  to  be  solved  consists  of  three/four 
equations:  continuity  and  two/ three  components  of  momentum  equation  in 
three/four  dependent  variables:  p,  u,  v,  w.  Using  notation  similar  to  that 
in  [15],  at  a  single  grid  point  this  system  of  equations  can  be  written  in 
the  following  form: 

3H((>)/3t  -  D((J)  +  S(*)  (1) 

where  ^  is  the  column-vector  of  dependent  variables,  H  and  S  are  column- 
vector  algebraic  functions  of  and  D  is  a  column  vector  whose  elements  are 
the  spatial  differential  operators  which  generate  all  spatial  derivatives 
appearing  in  the  governing  equation  associated  with  that  element. 

The  solution  procedure  is  based  on  the  following  two-level  implicit 
time-difference  approximations  of  (3): 

(Hn+l-  .  (d"*^)  +  (1-S)  (d"  +  S")  (2) 

where,  for  example,  denotes  and  At  ■  t"'*’^  -  t".  The 

parameter  8  (0.5  -  8  -  1)  permits  a  variable  time-centering  of  the  scheme, 
with  a  truncation  error  of  order  [At^,  (g  -  1/2)  At]. 

A  local  time  linearization  (Taylor  expansion  about  of  requisite 
formal  accuracy  is  introduced,  and  this  serves  to  define  a  linear  differen¬ 
tial  operator  L  (cf.  [15])  such  that 

-  d"  ♦  l”(*"*^-  •")  ♦  O(At^)  (3) 
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similarly 


rH+I  _  OH/a*)"  -  ♦")  +  0  (At^)  (4) 

gn+1  _  gti^  OS/a*)"  -  (j)">  +  0  (At^)  (5) 

Eqs.  (5-7)  are  inserted  into  Eq.  (4)  to  obtain  the  following  system  which  is 
linear  in 

(A  -  BAt  l")  -  4>")  -  At  (d"  +  s")  (6) 

and  which  is  termed  a  linearized  block  implilcit  (LBl)  scheme.  Here,  A 
denotes  a  matrix  defined  by 

A  i  (an/a*)"  -  BAt  (as/a(j.)*'  (7) 


Eq.  (8)  has  0  (At)  accuracy  unless  H  =  in  which  case  the  accuracy  is  the 
same  as  Eq.  (4). 

Special  Treatment  of  Diffusive  Terms 

The  time  differencing  of  diffusive  terms  is  modified  to  accomodate 
cross-derivative  terms  and  also  turbulent  viscosity  and  artificial  dissipa¬ 
tion  coefficients  which  depend  on  the  solution  variables.  Although  formal 
linearization  of  the  convection  and  pressure  gradient  terms  and  the  resulting 
implicit  coupling  of  variables  is  critical  to  the  stability  and  rapid  con¬ 
vergence  of  the  algorithm,  this  does  not  appear  to  be  important  for  the 
turbulent  viscosity  and  artificial  dissipation  coefficients.  Since  the 
relationship  between  and  dj  and  the  mean  flow  variables  is  not  conven¬ 
iently  linearized,  these  diffusive  coefficients  are  evaluated  explicitly  at 
t**  during  each  time  step.  Notationally,  this  is  equivalent  to  neglecting 
terms  proportional  to  3  or  3dj/3^  in  L**,  which  are  formally  pre¬ 

sent  in  the  Taylor  expansion  (5),  but  retaining  all  terms  proportional  to 
Vq  or  dj  in  both  L**  and  O'*. 

It  has  been  found  through  extensive  experience  that  this  has  little  if 
any  effect  on  the  performance  of  the  algorithm.  This  treatment  also  has  the 
added  benefit  that  the  turbulence  model  equations  can  be  decoupled  from  the 
system  of  mean  flow  equations  by  an  appropriate  matrix  partitioning  (cf. 

[IS])  and  solved  separately  in  each  step  of  the  ADI  solution  procedure.  This 
reduces  the  block  size  of  the  block  tridiagonal  systems  which  must  be  solved 
in  each  step  and  thus  reduces  the  computational  labor. 

In  addition,  the  viscous  terms  in  the  present  formulation  include  a  num 
cross-derivative  terms  implicitly  within  the  ADI  treatment  which  follows,  it 
is  not  at  all  convenient  to  do  so;  and  consequently,  all  cross-derivative 
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terms  are  evaluated  explicitly  at  t^.  For  a  scalar  model  equationrepresenting 
combined  convection  and  diffusion,  it  has  been  shown  by  Beam  and  Warming  that  the 
explicit  treatment  of  cross-derivative  terms  does  not  degrade  the  unconditional 
stability  of  the  present  algorithm.  To  preserve  notational  simplicity,  it  is 
understood  that  all  cross-derivative  terms  appearing  in  L**  are  neglected  but 
are  retained  in  O'*.  It  is  important  to  note  that  neglecting  terms  in  has 
no  effect  on  steady  solutions  of  Eq.  (8),  since  =  0  and  thus  Eq.  (8) 

reduces  to  the  steady  form  of  the  equations:  qh  4.  gn  _  q.  Aside  from 
stability  considerations,  the  only  effect  of  neglecting  terms  in  L"  is  to 
introduce  an  0  (At)  truncation  error. 

Consistent  Splitting  of  the  LBI  Scheme 

To  obtain  an  efficient  algorithm,  the  linearized  system  (8)  is  split 
using  ADI  techniques.  To  obtain  the  split  scheme,  the  multidimensional 
operator  L  is  revnritten  as  the  sum  of  three  "one-dimensional"  sub-operators 
(i  >  1,  2,  3)  each  of  which  contains  all  terms  having  derivatives  with 
respect  to  the  i-th  coordinate.  The  split  form  of  Eq.  (8)  can  be  derived 
either  as  in  [19,  22]  by  following  the  procedure  described  by  Douglas  and 
Gunn  [23]  in  their  generalization  and  unification  of  scalar  ADI  schemes,  or 
using  approximate  factorization.  For  the  present  system  of  equations,  the 
split  algorithm  is  given  by 

(A  -  BAtLj)  (<|»*  -  ♦")  -  At  (d"  +  S")  (8a) 

(A  -  BAtLj)  (♦**  -  ♦")  -  A  (♦*  -  ((>")  (8b) 

(A  -  8AtL")  (4»"*^  -  ♦”)  -  A  (♦**  -  ♦")  (8c) 

where  and  are  consistent  intermediate  solutions.  If  spatial 
derivatives  appearing  in  L£  and  D  are  replace  by  three-point  difference 
formulas,  as  indicated  previously,  then  each  step  in  Eqs.  (lOa-c)  can  be  solved 
by  a  block-tridiagonal  elimination. 

Combining  Eqs.  (lOa-c)  gives 

(A  -  6AtL“)  a"^  (a  -  BAtLp  a"^  (A  -  BAtl”)  -  ♦")  -  At  (d”  ♦  s”)  (9) 

which  approximates  the  tmsplit  scheme  (8)  to  0  (At^).  Since  the  intermediate 
steps  are  also  consistent  approximations  for  Eq.  (8),  physical  boundary 
conditions  can  be  used  for  ♦*  and  ♦**  [19,  22].  Finally,  since  the 
are  homogeneous  operators,  it  follows  from  Eqs.  (lOa-c)  that  steady  solutions 


51 


T 


have  the  property  that  ^  and  satisfy 

D“  ♦  S“  ■  0  (10) 

The  steady  solution  thus  depends  only  on  the  spatial  difference  approximations 
used  for  (12),  and  does  not  depend  on  the  solution  algorithm  itself. 
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SHOCK  LOCATION 


FIGURE  6-  CALCULATION  COORDINATE  SYSTEM  AND  CONSTANT  STATIC  PRESSURE 
CONTOURS  FOR  THE  CONVERGED  NORMAL  SHOCK  WAVE  /  BOUNDARY 
LAYER  INTERACTION  IN  A  TUBE 


FIGURE  7  “  PRESSURE  CONTOURS  CLOSE  TO  THE  SHOCK 


VELOCITY  VECTORS  IN  THE  NEAR  WALL  REGION  OF  THE  SHOCK  WAVE/ 
BOUNDARY  LAYER  INTERACTION 
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FIGURE  II-  PART  OF  THE  COORDINATE  SYSTEM  FOR  THE  LARGE  RADIUS 
OF  CURVATURE  BUMP  WITH  THIN  INLET  BOUNDARY  LAYER 
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FIGURE  16  -  CONSTANT  MACH  NUMBER  CONTOURS’  LARGE  RADIUS  OF 
CURVATURE  BUMP.  THICK  APPROACH  BOUNDARY  LAYER 


FIGURE  1 7  -  VELOCITY  VECTORS »  LARGE  RADIUS  OF  CURVATURE  BUMP, 
THICK  APPROACH  BOUNDARY  LAYER 


FIGURE  20-  VELOCITY  VECTORS*  SMALL  RADIUS  OF  CURVATURE  BUMP, 
THIN  APPROACH  BOUNDARY  LAYER 
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DOWNSTREAM  STATIC  PRESSURE  PERTURBATION 

PRESSURE  RESPONSE  AT  THE  ORIGINAL  SHOCK  LOCATION,  z/c  a  0.75 
STATIC  PRESSURE  RESPONSE  AT  z/c  =  0.61 
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FIGURE  22  -  TIME  HISTORY  OF  WALL  STATIC  PRESSURE  FOR  TRANSONIC,  TURBULENT 

FLOW  OVER  AN  AXISYMMETRIC  BUMP  AT  REDUCED  FREQUENCY  (f^  =  0.175) 


FREQUENCY  RESPONCE  FOR  UNSTEADY  FLOW  OVER  AN  AXISYMMETRIC 
BUMP  {f.  =  0.004). 


75 


FIGURE  25-  TIME  HISTORY  OF  CONSTANT  STATIC  PRESSURE  CONTOURS  FOR  AXISYMMETRIC 
BUMP  FLOW  SUBJECT  TO  OSCILLATING  DOWNSTREAM  PRESSURE 
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FIGURE  26—  COMPRESSION  RAMP-  COORDINATE  SYSTEM 


FIGURE  27  -  COMPRESSION  RAMP-  DENSITY  CONTOURS 


FIGURE  28  -  COMPRESSION  RAMP=  PRESSURE  CONTOURS 
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FIGURE  29  -  STATIC  PRESSURE  DISTRIBUTION*  AXISYMMETRIC  COMPRESSION  RAMP 


FIGURE  31-  CONSTANT  PRESSURE  CONTOURS  FOR  SHOCK  IMPINGING  ON 
A  TURBULENT  BOUNDARY  LAYER,  USING  SHOCK  ADAPTING 
COORDINATE  SYSTEM,  M  =  2.96,  SHOCK  GENERATOR  ANGLE 
10.89  DEGREES. 
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